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Abstract 

This  document,  together  with  the  software  and  test  datasets  supplied 
and  installed  on  a  SUN  workstation  at  Phillips  Laboratory  are  the  com¬ 
pletion  of  Task  1  of  the  contract  3-D  PIC  Software  for  MIMD  Computers, 
PIIN;F61708-93-C0011. 

1  INTRODUCTION 

The  contents  of  this  document  and  its  Appendices  and  Annexes  describe  the 
software  and  test  datasets  for  Task  1  of  the  contract  S-D  PIC  Software  for 
MIMD  Coniputers,  PIIN:F61708-93-C0011.  Further  details  on  the  algorithms 
and  software  design  are  given  in  earlier  documents  issued  during  this  contract 
[10,  11,  12,  13,  14,  16],  and  during  earlier  studies[3,  4], 

The  objectives  of  the  work,  as  stated  in  the  contract,  are 

Task  1:  The  primary  target  for  the  present  proposal  is  to  implement  a  fully 
three  dimensional  time  domain  electromagnetic  field  solver  code  using 
body-fitted  brick  finite  elements.  The  code  will  be  Fortran  based,  and 
designed  to  run  efficiently  on  distributed  memory  MIMD  computers. 

Task  2:  Task  2  will  extend  the  code  developed  in  Task  1,  adding  particle  in¬ 
tegration  modules  to  the  electromagnetic  solver  to  create  a  3-D  general 
geometry  PIC  code. 

As  described  in  the  original  proposal  [5],  the  method  of  achieving  these  ob¬ 
jectives  has  been  based  on  that  used  for  the  two  dimensional  benchmarking 
program  [4]  developed  under  Air  Force  Office  of  Scientific  Research  (AFSC) 
Contract  F49620-92-C-0035.  Only  Task  1  was  scheduled  to  be  undertaken  in 
the  present  phase  of  the  work,  although  some  of  the  Task  2  items  have  been 
started  to  help  achieve  the  goals  of  Task  1. 

The  numerical  algorithms  used  in  the  electromagnetic  software  are  based 
on  those  described  in  Reference  [3].  In  summary,  these  use 
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1.  the  ‘Virtual  Particle’  derivation  method  [2]  applied  to  tensor  field  com¬ 
ponents, 

2.  a  multiblock  spatial  decomposition  applied  to  both  fields  and  particles, 

3.  transfinite  interpolation  subdivision  of  the  curvilinear  hexahedral  multi¬ 
blocks  into  hexahedral  elements, 

4.  indirect  (“glue  patch”)  addressing  between  multiblocks,  and  logical  cubi¬ 
cal  mesh  node  addressing  within  blocks, 

5.  lumped  approximations  to  the  finite  element  equation  [3]. 

The  initially  proposed  scheme  has  evolved  a  little  in  the  light  of  experience,  as 
will  be  apparent  from  the  material  presented  in  this  report  and  its  Annexes. 

Task  1:  Electromagnetic  Software  The  stated  goals  of  the  first  part  of 
the  contract  were  [5]  : 

Subtask  1.1  Identify  and  evaluate  the  mesh  generation  software  to  be  used. 
Specify  its  interface  to  the  electromagnetic  software. 

Subtask  1.2  Prepare  a  document  containing  the  specification  of  the  software 
and  validation  test  cases.  Agree  specifications  with  USAF  technical  con¬ 
tract  monitor. 

Subtask  1.3  Write  the  3-D  electromagnetic  solver  program. 

Subtask  1.4  Execute  test  cases  on  SUN  SPARCstation  and  on  the  agreed 
distributed  memory  MIMD  computer. 

Subtask  1.5  Deliver  documented  software  and  test  run  input  and  output  to 
Phillips  laboratory. 

Subtasks  1.1  and  1.2  have  been  documented  in  earlier  reports  [10.  11,  12,  13] 
issued  under  this  contract.  Annexes  B,  C,  D,  F  and  G  update  and  extend  the 
specifications.  The  remaining  activities  undertaken  to  satisfy  the  requirements 
of  Subtasks  1..3-1.5,  and  to  undertake  some  preparatory  work  for  Task  2  are 
outlined  in  this  document  and  its  Annexes. 

The  .Appendices  to  this  report  contain  a  list  of  the  program  units,  and 
some  fragments  of  the  units  and  integration  test  outputs  used  in  developing 
and  verifying  the  code.  The  Annexes  contain  further  analysis  and  specification 
documents,  plus  copies  of  two  papers  arising  from  the  work  undertaken.  These 
Annexes  are: 

Annex  A: 

R  W  Hockney, 

LPM3  Benchmark  ResiiHs  on  the  Intel  Paragon  and  iPSC/860, 

Technical  Note  AEA/TYKB/31878/TN/6, 

Culham  Laboratory,  July  1994. 
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A  principal  objective  of  the  work  has  been  to  implement  the  body  fitted 
PIC  software  on  parallel  computers.  The  benchmark  code  has  provided 
the  testbed  for  the  development  of  the  interblock  message  passing  pro¬ 
gram  units  for  the  Paragon  and  iPSC.  This  Note  summarises  the  parallel 
performance  obtained  for  a  simple  3-D  electron  plasma  configuration. 

Annex  B: 

W  Arter, 

The  System  of  Dimensionless  Units, 

Technical  Note  AEA/TYKB/31878/TN/7, 

Culham  Laboratory,  September  1994. 

The  originally  proposed  system  of  internal  units  described  in  Ref  [3]  has 
been  revised.  This  Note  outlines  the  new  choice  of  units  for  the  kernel  of 
the  simulation  program  PIC3D. 

Annex  C: 

W  Arter, 

Dispersion  and  Stability  Analysis  for  Maxwell’s  Equations, 

Technical  Note  AEA/TYKB/31878/TN/8, 

Culham  Laboratory,  September  1994. 

This  Note  derives  the  stability  criterion  used  to  specify  the  timestep  used 
for  integration  of  the  field  equations  in  PIC3D. 

Annex  D: 

N  J  Brealey,  J  W  Eastwood  and  W  Arter, 

SDPIC  Diagnostics, 

Technical  Note  AEA/TYKB/31878/TN/9, 

Culham  Laboratory,  September  1994. 

This  Note  presents  the  diagnostics  u.sed  by  PIC3D. 

Annex  E: 

N  J  Brealey, 

MPICTIM  User’s  Guide, 

Technical  Note  AEA/TYKB/31878/TN/10, 

Culham  Laboratory,  September  1994. 

The  timeseries  output  files  (.tsd  files)  produced  by  PIC3D  diagnostics 
can  be  examined  on  a  UNIX  Workstation  using  the  OSF/Motif  GUI  tool 
MPICTIM.  This  Note  describes  how  to  use  MPICTIM. 

Annex  F: 

W  Arter, 

Boundary  Conditions  for  Maxwell’s  Equations  in  General  Geometry, 
Technical  Note  AEA/TYKB/31878/TN/11, 

Culham  Laboratory,  September  1994. 

This  Note  contains  a  derivation  of  boundary  conditions  in  general  geom¬ 
etry,  and  a  discussion  of  their  implementation  in  PIC3D. 
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Annex  G: 

W  Arter, 

Extended  Description  of  Inputs  used  to  Generate  Meshes, 

Technical  Note  AEA/TYKB/31878/TN/13, 

Culhani  Laboratory,  September  1994. 

This  Note  extends  the  specification  give  in  Ref  [13],  and  supersedes  that 
document. 

Annex  H: 

J  W  Eastwood,  W  Arter,  N  J  Brealey  and  R  W  Hockney, 

Body  Fitted  PIC  Software  for  Microwave  Device  Modelling, 

EUROEM,  Bordeaux,  France  May  1991. 

Annex  I: 

\V  Arter,  J  \V  Eastwood,  N  J  Brealey  and  R  W  Hockney, 
Electromagnetic  Modelling  in  Arbitrary  Geometry  by  PIC  Methods  on 
MIMD  Computers, 

pp  297-300  in  6th  Joint  EPS-.4PS  Int  Conf  on  Phys  Computing  PC’94 
(R  Gruber  and  M  Tomassini,  eds),  EPS,  Geneva(1994). 

The  original  proposal  was  based  on  a  bottom-up  development  plan,  where 
the  aim  was  to  develop  a  working  kernel  where  the  geometry  could  only  be 
changed  by  some  recoding,  and  then  to  add  the  user  functionality  later.  The 
design  stage  of  the  programme  of  work  for  Task  1  led  to  a  substantial  revision 
of  this  plan;  the  decision  w’as  made  to  reschedule  the  software  development  into 
a  top-downwards  scheme,  with  unit  testing  of  groups  of  program  units. 

The  top-downwards  approach  begins  with  the  preprocessor’s  user  input  in¬ 
terface.  then  develops  code  to  validate  input,  to  display  the  data,  and  to  trans¬ 
form  the  data  into  the  form  needed  by  the  simulation  kernel.  This  pre])rocessor 
and  its  interface  to  the  kernel  have  proved  far  more  difficult  and  time  consum¬ 
ing  than  envisaged.  The  initial  plan  to  use  the  EAGLE  grid  generation  package 
[8]  were  abandoned  after  evaluating  E.VGLE;  our  assessment  of  that  software 
indicated  that  it  could  be  of  great  value  in  our  electromagnetic  code,  but  to 
exploit  it  would  require  an  investment  of  effort  that  was  beyond  the  scope  of  the 
present  contract.  Consequently,  we  have  replaced  EAGLE  input  by  a  less  ambi¬ 
tious  grid  generation  capability,  but  have  left  open  the  possibility  of  importing 
element  nets  from  it  at  a  later  date. 

The  need  for  data  validation  and  graphical  display  led  to  a  number  of  mod¬ 
ules  being  moved  from  the  simulation  kernel  to  the  preprocessor,  which  in  turn 
has  required  substantial  reworking  of  the  data  interface.  In  addition,  the  delay 
in  getting  adequate  data  initialisation  has  hampered  the  development  and  de¬ 
bugging  of  the  electromagnetic  routines  of  the  main  simulation  program.  The 
design  is  now  stable,  but  not  complete  in  all  its  details. 

The  project  has  been  a  major  software  development  effort,  with  more  than 
one  new  subprogram  being  created  per  day,  and  with  the  same  number  of 
existing  program  units  being  adapted  and  integrated  into  the  developing  code 
structure.  The  subprograms  have  been  organised  by  function  into  a  number 
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of  libraries,  each  under  SCCS  configuration  control  to  assure  quality.  These 
libraries  are  described  below  in  Section  3  and  in  Appendix  A. 

2  SOFTWARE  DEVELOPMENT 

A  history  of  the  development  of  the  software  is  given  by  the  monthly  Progress 
Statements  [15]  issued  to  the  Contract  Monitor.  Inspection  of  those  documents 
reveals  a  four  part  parallel  development:  preprocessing,  simulation,  diagnostics 
and  distributed  memory  MIMD  computer  implementation.  The  first  three  of 
these  are  closely  coupled  by  the  need  to  develop  compatible  interfaces.  The  final 
item,  the  development  of  the  message  passing  for  parallel  implementation,  has 
largely  been  isolated  by  using  LPM3,  a  3-D  extension  of  the  2-D  benchmark¬ 
ing  program  LPM2.  LPhI3  has  provided  a  test-bed  for  new  message  passing 
software  for  use  on  the  iPSC  and  Paragon  computers  (see  Annexes  A  and  H). 
The  four  parts  of  the  software  development  correspond  to  the  four  codes 


Code 

Purpose 

PEGGIE 

PIC3D 

MPICTIM 

MIMDPIC(LPM3) 

preprocessing 

simulation 

time  series  data  analysis 

3-D  parallel  benchmarking 

MPICTIM  is  an  interactive  program  which  uses  a  Motif  GUI  interface.  It  is 
a  general  purpose  time  series  data  analysis  program  which  was  not  developed 
under  this  contract,  and  is  only  being  provided  in  an  executable  form  for  use 
on  SUN  workstations.  The  large  number  of  program  units  shared  by  the  other- 
codes  are  split  into  11  libraries,  organised  by  function: 


Library 

Purpose 

OLYMPUS 

CRONUS 

NETBLK 

NETGLB 

EMBLK 

EMGLB 

PARBLK 

PARGLB 

GLBLIB 

MCDLIB 

DIAG 

OLYMPUS  utility  library 

OLYMPUS  template  library 

block  element  net  library 

global  element  net  library 

block  electromagnetic  library 

global  electromagnetic  library 

block  particle  library 

global  particle  library 

global  data,  and  related  program  units 

machine  dependent  library 

diagnostic  output  library 

The  subprograms  included  in  the  codes  and  libraries  and  a  brief  description  of 
their  functions  are  listed  in  Appendix  A. 

The  main  effort  has  been  directed  towards  implementing  those  parts  of 
the  software  required  to  initialise  and  execute  the  electromagnetic  part  of  the 
calculation,  although  some  effort  (particularly  in  the  benchmarking  code)  has 
been  devoted  to  implenrenting  the  core  of  the  particle  integration  software. 

The  following  Sections  give  brief  descriptions  of  the  codes  and  code  modules. 
Section  5  summarises  the  test  datasets  provided  with  the  software. 
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3  DESCRIPTION  OF  CODES 

The  suite  of  codes  stored  in  directory  3DPIC  consists  of  the  preprocessor,  the 
main  simulation  code,  the  postprocessor  and  the  parallel  benchmarking  code. 

3.1  Preprocessing  (PEGGIE) 

The  concept  underlying  the  PEGGIE  (Parallel  Electromagnetic  General  Geom¬ 
etry  Interface)  approach  to  grid  generation,  boundary  conditions  and  material 
property  specification  is  that  the  device  to  be  modelled  is  made  by  joining 
‘parts’  together.  An  important  advantage  is  the  a  ‘part’  might  be  output  of  an¬ 
other  geometry  package,  although  normally  it  will  consist  of  a  subroutine  and 
the  associated  variables.  A  simple  example  of  the  specification  of  parts  is  given 
in  Appendix  B;  there  a  section  of  rectangular  waveguide  is  specified  in  terms 
of  two  blocks,  placed  side  by  side  in  the  2-direction  and  with  conducting  walls 
in  the  1-  and  .3-directions,  and  periodic  conditions  in  the  2-direction.  Further 
details  are  given  by  the  specification  provide  in  Annex  G. 

Input  to  PEGGIE  consists  of  tagged  records  each  followed  by  FORTRAN 
NAMELIST  input  to  assign  parameters  within  the  record,  which  may,  for  ex¬ 
ample.  define  the  geometry  of  a  part  of  block  of  a  device.  Within  the  limitations 
of  the  command  line  format  PEGGIE  is  designed  to  be  user  friendly,  ie  it  tol¬ 
erates  certain  mistakes,  tries  to  fill  in  gaps  in  input  data  and  will  flag  some 
unusual,  but  not  necessarily  incorrect,  states.  There  is  extensive  validation 
of  input  parameters,  wherein  grid  spacings  may  be  adjusted  to  make  different 
parts  fit. 

The  second  phase  of  PEGGIE  is  concerned  with  generating  the  surface  and 
line  patches  that  describe  block  connectivity.  Global  co-ordinates  for  the  part 
boundaries  are  also  produced.  Copious  diagnostics  for  the  developer’s  benefit 
are  also  generated. 

For  other  users  the  main  output  consists  of  projections  of  the  proposed 
device  (to  check  that  it  has  been  specified  correctly)  and  a  data-file  that  serves 
as  input  to  the  PIC3D  code.  Thus  PEGGIE  is  a  general  interface  in  the  sense 
that  control  parameters  and  initial  conditions  for  PIC3D  may  be  specified. 

PEGGIE  uses  routines  from  the  NETGLB,  NETBLK,  EMBLK,  GLBLIB, 
CRONUS  and  OLYMPUS  libraries.  There  are  currently  218  program  units, 
including  those  from  the  shared  libraries.  Figure  1  shows  the  interdependence 
of  the  libraries.  The  links  to  OLYMPUS  are  omitted  for  clarity;  all  other 
program  modules  use  the  OLYMPUS  module  library. 

3.2  Simulation  (PIC3D) 

The  simulation  code,  PIC3D,  is  kept  as  small,  simple  and  as  fast  as  possible. 
Nearly  all  the  data  initialisation  is  performed  by  the  preprocessor  PEGGIE, 
and  as  much  as  possible  of  the  data  processing  for  diagnostics  is  performed 
by  postprocessors.  Each  process  in  the  parallel  implementation  of  the  main 
code  writes  separate  output  files  which  for  simplicity  are  concatenated  during 
postprocessing. 
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Figure  1:  A  schematic  of  the  library  dependencies  for  the  preprocessor  program 
PEGGIE 


The  main  code  uses  dynamic  storage  allocation  so  that  the  code  does  not 
need  to  be  recompiled  for  different  sized  problems  and  can  use  the  minimum 
amount  of  storage  required  for  each  problem.  Dynamic  load  balancing  has  not 
been  implemented,  but  could  be  straightforwardly  added  without  reworking  the 
data  structures  used. 

All  the  data  control  logic  and  a  number  of  particle  integration  subprograms 
are  already  included  in  PIC3D,  although  at  present  these  routines  will  work 
only  for  simple  geometries.  These  will  be  developed  further  under  the  Task  2 
programme  of  work. 

The  main  timestep  loop  has  been  structured  to  include  electromagnetic 
field  integration  subcycling  with  respect  to  the  particle  integration.  Different 
timestepping  for  different  blocks,  although  possible,  has  not  been  implemented. 
However,  some  relaxation  of  timestep  restriction  in  certain  cases  may  be  possible 
by  using  the  FFT  filtering  routines  which  have  been  incorporated  into  the  EM 
solver  for  use  with  cylindrical  blocks;  these  routines  are  intended  for  use  both 
as  noise  filters  and  means  of  relaxing  the  Courant  condition  near  cylindrical 
mesh  singularities. 

PIC3D  uses  routines  from  the  NETGLB,  NETBLK,  PARGLB,  PARBLK, 
EMGLB,  EMBLK,  GLBLIB,  MCDLIB,  CRONUS  and  OLYMPUS  modules. 
It  currently  uses  about  100  program  units.  Figure  2  shows  the  main  program 
module  and  the  interdependence  of  the  support  libraries.  The  links  to  OLYM¬ 
PUS  are  omitted  for  clarity:  all  program  modules  depend  on  the  OLYMPUS 
library. 
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(^NETBLK^  EMBLK  )  C  GLBLIB  )  (  paRblk) 


Figure  2:  The  main  simulation  code  PIC3D  and  its  support  libraries.  The 
shaded  libraries  are  those  being  developed  during  the  Task  1  phase  of  the 
project. 


The  PIC3D  code  could  generate  an  enormous  amount  of  graphical  data. 
When  analysing  a  simulation  it  is  vital  that  data  from  the  simulation  is  filtered 
to  reduce  the  quantity  of  data,  and  data  is  stored  in  the  most  compact  form 
possible. 

Two  classes  of  snapshot  output  are  possible: 

1.  Line  plots  show  the  dependence  of  a  scalar  quantity  on  the  axial  position 
along  the  device.  The  plots  are  usually  simple  x-y  line  graphs. 

2.  Slice  plots  show  the  dependence  of  a  scalar  quantity  on  a  slice  through 
the  device.  The  plots  are  usually  contour  plots  or  false  colour  plots. 

The  data  input  specification  to  produce  these  snapshot  outputs  and  to  produce 
datasets  for  MPICTIM,  together  with  a  summary  of  the  possible  outputs,  is 
given  in  Annex  D. 

The  main  simulation  code  takes  data  for  selected  plots  and  saves  the  data 
in  tagged  unformatted  files  for  subsequent  analysis.  Diagnostics  in  the  main 
code  may  integrate  over  surfaces,  average  over  time  intervals  or  perform  other 
processing  before  the  data  is  saved.  The  snapshot  data  includes  data  which  is 
averaged  over  a  time  interval  as  well  as  instantaneous  values. 

Particle  plotting  routines  presently  included  in  the  libraries  only  work  for 
simple  cartesian  test  cases.  The  will  be  replaced  by  more  robust  general  geom¬ 
etry  routines  in  the  course  of  the  Task  2  programme  of  work. 

3.3  Timeseries  Analysis  (MPICTIM) 

Another  output  from  the  PIC3D  is  timeseries  data.  This  data  is  stored  in 
files  using  the  same  format  as  the  2D  code  TWT  and  2-D  MILO  simulation 
codes  developed  at  Culham.  although  the  data  items  stored  in  the  file  may  be 
different.  The  timeseries  analysis  tools  may  be  used  to  analyse  data  from  the 
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3-D  code.  In  particular,  MPICTIM,  the  timeseries  analysis  tool  with  a  Motif 
GUI,  may  be  used. 

A  user’s  guide  for  MPICTIM  is  contained  in  Annex  E.  MPICTIM  was  not 
developed  under  this  contract,  but  an  executable  for  the  SUN  workstation  is 
provided  to  view  the  timeseries  data  generated  by  PIC3D  test  runs. 

3.4  Parallel  Benchmarking  Code  (MIMDPIC/LPM3) 

LPM3  is  a  version  of  the  main  simulation  code  PIC3D  with  built  in  initiali¬ 
sation  routines,  particle  integration  routines,  restricted  geometry,  and  limited 
output.  It  can  perform  the  test  cases  summarised  in  [4],  and  model  3-D  elec¬ 
tron  plasmas  in  rectangular  brick  regions.  It  was  developed  from  the  earlier 
2-D  benchmarking  code  [4].  The  serial  Implementation  in  directory  MIMDPIC 
is  used  to  generate  test  datasets  for  PIC3D  and  to  provide  cross  checking  data 
for  the  parallel  implementation  running  on  the  Sandia  Paragon  machine. 

The  version  of  LPM3  on  the  Workstation  uses  the  various  support  libraries 
for  many  of  its  program  units.  The  LPM3  benchmark  version  installed  on  the 
Sandia  Paragon  XP/S  140  under  the  directory  "rhockne/Paragon/LPM3  was 
used  to  obtain  the  results  from  LPM3  which  are  presented  in  Annexes  A,  II 
and  I. 

4  DESCRIPTION  OF  LIBRARIES 

This  Section  gives  a  description  of  each  of  the  support  library  modules  used 
by  the  suite  of  codes  for  3-D  modelling.  There  are  currently  11  modules  which 
are  used  by  the  preprocessor  and  the  main  code.  There  are  a  large  number  of 
program  units,  test  units,  COMMON  and  documentation  files.  The  routines 
in  each  module  are  listed  in  Appendix  A;  further  documentation  is  included  in 
the  .doc  files  in  each  of  the  program  and  library  subdirectories. 

The  OLYhIPUS  and  CRONUS  libraries  contain  routines  used  to  provide  a 
standard  software  development  environment  [1,  chap  4].  The  remaining  libraries 
are  specific  to  the  electromagnetic  Particle-in-Cell  software. 

The  three  main  parts  of  the  calculation  are  the  element  net  generation 
(NET...),  the  electromagnetic  field  calculation  (EM...)  and  the  particle  inte¬ 
gration  (PAR... ).  Program  units  in  each  of  the  categories  are  grouped  into  block 
(...BLK)  and  global  (...GLB)  subprograms. 

Block  subprograms  operate  on  data  within  a  single  curvilinear  hexahedral 
block,  and  do  not  have  access  to  the  global  connectivity  data,  of  the  multiblock 
decomposition.  NETBLK,  EMBLK  and  PARBLK  are  block  libraries. 

Global  subprograms,  which  handle  the  interconnection  of  the  mulitblocks, 
use  the  block  subprograms  as  primitives  for  data  manipulation.  NETGLB, 
EMGLB  and  PARGLB  are  global  libraries.  The  global  net  data  manipulated 
by  these  routines  is  either  held  in  the  COMMON  blocks  defined  in  module 
GLBLIB,  or  in  the  case  of  the  principal  field,  particle  and  buffering  arrays,  in 
dynamically  allocate  memory. 

The  mapping  of  global  data  onto  processes  and  from  processes  onto  proces¬ 
sors  is  machine  dependent.  Routines  for  these  operations,  which  will  differ  on 
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different  architecture  machines,  are  collected  into  the  MCDLIB  library.  Each 
subdirectory  of  the  3DPIC  directory  (apart  from  MPICTIM)  contains  a  library 
or  program  also  contains  a  number  of  documentation  modules,  FORTRAN 
source  and,  in  some  cases  common  blocks. 

4.1  CRONUS 

This  library  contains  the  basic  template  for  OLYMPUS  codes  [1,  9].  The  library 
contains  dummy  versions  of  routines  which  are  normally  replaced  by  routines 
with  the  same  name  in  the  main  module.  These  routines  were  previously  in  the 
OLYMPUS  library  in  earlier  versions  of  the  OLYMPUS  environment. 

The  subroutine  BASIC  has  been  enhanced  so  that  codes  fail  gracefully  if 
they  can  not  open  the  files  they  need. 

4.2  OLYMPUS 

This  library  contains  the  OLYMPUS  utility  routines  (CYCLOPS)  [9].  New 
utility  routines  have  been  added  for  reading  the  command  line  (CMDLIN), 
concatenation  of  blank  padded  character  strings  (CONCAT).  The  *VAR  and 
*ARRAY  routines  have  been  updated  to  allow  variable  length  labels  and  new 
*OUT  routines  have  been  added  for  writing  datasets.  The  UN  AMES  routine 
provides  system  information. 

4.3  NETBLK 

This  module  contains  routines  for  element  NETs  in  a  BLocK.  It  contains  rou¬ 
tines  for  generating  meshes  and  routines  u.sed  to  visualise  meshes.  The  test  unit 
TEST  tests  various  different  types  of  block  geometries. 

4.4  NETGLB 

This  module  contains  routines  for  global  net  assembly.  The  test  unit  TEST 
tests  connecting  different  types  of  blocks  together. 

4.5  EMBLK 

This  library  contains  the  routines  for  the  ElectroMagnetic  solver  in  a  BLocK. 
There  are  three  different  versions  of  the  Ampere  ecpiation  routine: 

1.  AMPERE  is  a  simple  implementation  which  requires  separate  storage  for 
H,  j  and  d.  It  loops  through  components  of  d  outside  the  loop  which  loops 
through  nodes.  It  uses  the  standard  one  dimensional  array  addressing  [3] 
which  allows  the  same  compact  code  to  do  2D  and  3D  calculations. 

2.  AMPINM  is  an  implementation  which  overwrites  H  with  d,  uses  mul¬ 
tidimensional  arrays,  and  requires  two  planes  of  bordering  elements  as 
workspace. 
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3.  AMPINV  is  the  same  as  AMPINM  except  that  it  is  for  the  outer  regions 
where  there  is  no  particle  current.  In  the  absence  of  particle  or  other 
volume  filling  currents,  it  allows  reduced  storage  by  eliminating  the  need 
for  current  array  space  to  be  reserved  for  the  block. 

The  routines  GBTOH  and  GDTOE  convert  b  to  H  and  d  to  E  respectively. 
The  routine  GORTHV  sets  vacuum  block  metrics.  The  routine  GLUEIO  trans¬ 
fers  data  to  and  from  the  glue  patch  buffers.  Routines  for  embedded  conductors 
and  electromagnetic  boundary  conditions  are  contained  in  this  module.  These 
are  at  present  only  partially  developed. 

The  Fast  Fourier  Transform  routines  which  are  used  for  spectral  filtering 
are  in  this  module. 

4.6  EMGLB 

This  modules  contains  the  ElectroMagnetic  GLoBal  routines.  In  general,  these 
global  routines  access  global  data  in  COMMON  whereas  the  routines  in  EM- 
BLK  access  data  via  argument  lists. 

GETMET  generates  metrics,  BLKIO  copies  data  between  blocks  and  patches 
and  other  routines  deal  with  EM  boundary  conditions. 

4.7  GLBLIB 

This  is  the  GLoBal  data  storage  LIBrary.  The  main  units  in  this  module  are 
the  COMMON  blocks  for  data  storage.  It  also  contains  routines  which  initialise 
the  COMMON  blocks  and  read  or  write  their  contents  from  or  to  files. 

4.8  MCDLIB 

All  routines  (apart  from  those  in  the  OLYMPUS  library)  which  are  MaChine 
Dependent  are  grouped  into  this  library.  At  present,  these  are  primarily  the 
routines  which  exchange  information  between  patch  buffers.  Special  versions 
of  these  routines  implement  the  message  passing  in  the  benchmarking  program 
on  the  iPSC  and  Paragon  computers.  For  use  on  workstations,  the  message 
passing  simply  becomes  memory  to  memory  copying. 

4.9  PARBLK 

This  module  contains  the  local  particle  routines.  It  includes  routines  for  particle 
boundary  conditions.  The  routine  BEAMIN  is  used  for  injecting  a  prescribed 
beam.  It  uses  a  table  to  specify  the  beam  parameters  and  places  particles  on  the 
block  surface  and  gives  them  a  launch  time.  Support  routines  use  the  launch 
time  to  compute  the  particle  position  at  the  end  of  the  time  step.  The  routines 
for  injected  saved  particles  share  some  of  the  same  support  routines. 

At  present,  these  routines  have  been  developed  only  partially  for  the  pur¬ 
poses  of  the  benchmarking  tests  in  the  iPSC  and  Paragon  computers. 
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4.10  PARGLB 

This  module  contains  the  global  particle  routines.  It  includes  routines  for  load¬ 
ing  and  unloading  particle  buffers. 

4.11  DIAG 

This  module  contains  the  routines  which  collect  and  output  the  line,  line  inte¬ 
gral,  surface  and  volume  integral  diagnostics  as  described  in  Annex  D. 

5  TEST  DATA 

The  four  programs:  PEGGIE,  PIC3D,  MIMDPIC  and  MPICTIM  and  the 
eleven  associated  libraries  have  been  subjected  to  a  number  of  unit  test  and 
integration  tests  to  verify  the  coding.  The  ultimate  goal  of  the  development  is 
to  use  PEGGIE  (and  perhaps  other  data  preparation  tools)  to  handle  the  input 
geometry,  initial  values,  boundary  values,  material  properties,  numerical  and 
computation  parameters  and  output  selections.  PEGGIE  validates  this  data 
(the  User  Input  File  or  .uif  file)  and  prepares  the  simulation  input  datasets 
( .dat  file)  for  use  by  the  simulator  PIC3D.  PIC3D  outputs  a  number  of  snap¬ 
shot  files  and  timeseries  data  (.tsd)  files  for  analysis  by  MPICTIM  and  the 
data  analysis  and  visualisation  tools. 

Described  below  are  some  of  the  unit  tests  and  integration  tests  provided 
with  the  software  delivered  to  the  Phillips  Laboratory.  Datasets  appropriate 
for  testing  the  source  on  workstations  are  installed  in  the  appropriate  TEST 
subdirectories  under  the  directories  containing  the  library  or  main  code  sources. 

5.1  CRONUS 

The  unit  tests  of  CRONUS  (Appendix  D)  runs  the  default  OLYMPUS  skeleton 
program  and  tests  a  number  of  the  OLYMPUS  library  routines. 

5.2  OLYMPUS 

The  OLYMPUS  library  provides  support  for  all  the  programs,  and  so  is  tested 
in  all  the  unit  and  integration  tests.  A  unit  test  for  the  machine  dependent 
routine  is  included  in  the  library  and  is  involved  by  typing  eg 

make  test . solaris2 

for  testing  the  SUN  Solaris  2  implementation.  Reference  test  outputs  for  SUN 
Solaris  1  (see  Appendix  E),  SUN  Solaris  2,  Silicon  Graphic  Iris  and  IIP  720 
\^Urkstations  are  provided. 

5.3  NETBLK 

A  unit  test  program  xtest  is  contained  in  subdirectory  TEST  of  directory 
NETBLK.  To  execute  the  test  program,  typo 
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xtest  test.dat 

An  example  user-program  dialogue  is  shown  in  the  listing  of  the  test .  log  in 
Appendix  F.  In  the  test,  VPNX,  VPNY,  VPNZ  are  the  viewing  position  for  plotting 
the  block,  and  PLOTYPE  =  1,  2  and  3  respectively  select  the  plotting  of  all 
element  edges,  element  edges  on  the  block  surface  and  block  edges  only. 

The  input  dataset  test.dat,  the  output  file  o_test,  and  a  collage  of  the 
eight  block  plots  generated  (file  g_test)  are  also  included  in  Appendix  F. 

5.4  NETGLB 

Appendix  G  contains  a  unit  test  which  demonstrates  the  successful  assembly  of 
a  number  of  blocks  into  a  cylindrical  mesh.  This  test  demonstrates  the  correct 
logical  and  physical  net  assembly,  and  the  mesh  projection  output  routines. 
The  files  listed  are 

mtest.log  :  a  record  of  the  interactive  dialogue 
mtest  .dat  :  the  input  dataset 
o_mtest  :  printed  output 

g.mtest  :  graphical  output 

5.5  EMBLK 

The  principal  testing  of  the  electromagnetic  routines  is  undertaken  in  the  frame¬ 
work  of  the  programs  MIMPDIC  and  or  PIC3D.  However  there  are  four  unit 
tests  included  in  the  EMBLK  directory: 

1.  TEST  tests  the  FFT  routines 

2.  TESTl  tests  AMPINM  and  its  support  routines 

3.  TEST2  tests  AMPINV  and  its  support  routines 

4.  TEST3  will  test  the  field  solver  using  cylindrical  cavity  modes. 

Listings  of  output  from  these  tests  are  included  in  Appendix  II. 

5.6  PEGGIE 

There  are  3  test  cases  for  the  PEGGIE  code  in  the  TEST  directory.  Case  xxx 
runs  with  data  from  the  input  file  xxx.uif  and  generates  files  o_xxx,  xxx.dat 
and  g_xxx.  The  cases  are  designed  primarily  to  check  out  the  geometry  and 
normally  will  require  some  modifications  to  alter  output  and  other  parameters 
before  they  can  be  used  as  input  to  PIC3D.  In  particular,  mug2.dat  cannot 
sensibly  be  used  with  PIC3D,  since  the  handle  of  the  mug  is  only  one  element 
wide.  Case  tl30  models  a  cylindrical  cavity  made  by  assembling  orthogonal 
annular  and  polar  blocks.  Case  uni20  generates  a  rectangnlar  cavity  from  8 
non-orthogonal  blocks.  The  .uif  files  are  listed  in  Appendix  B. 

The  o_xxx  files  are  of  interest  solely  to  the  developer,  so  only  the  xxx.dat 
files  need  checking.  It  is  however  helpful  to  record  the  graphics  parameters 
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required  to  generate  the  g_xxx  files,  as  the  reproduction  of  these  is  a  simple 
and  powerful  test  of  code  validity.  Usually  only  the  plotting  window  needs 
changing  from  its  default  values.  For  the  cases  inug2,  tl30  and  uni20  the 
values  used  are  respectively  (-25,  5,  -15  -15),  (-25,  55,-40,  40)  and  (-10,  90,  -30, 
70).  For  inug2,  the  first  parameters  have  also  to  be  changed  to  (-0.3,  0.2,  0.7,  0, 
0,  0,-1,  0,0). 


5.7  MIMDPIC 

The  test  data  sets  provided  with  the  earlier  benchmarking  program  [4]  have 
been  updated  to  provide  test  data  for  the  evolving  3DPIC  codes.  Listed  below 
are  test  datasets  for  the  current  MIMDPIC  program.  Note  that  the  names 
have  been  changed  to  user  input  (.uif)  files,  although  there  is  a  one  to  one 
correspondence  with  the  .dat  files  in  Ref  [4]. 


testl .uif 
test2.uif 
tests .uif 
testd.uif 
tests .uif 
teste .uif 
test? .uif 
tests .uif 
tests .uif 
testlO.uif 
testll .uif 
testl2 .uif 
testis .uif 
testl4.uif 
testis .uif 
testie .uif 
testl7 .uif 
testis .uif 
testis .uif 
test20 .uif 
test21 .uif 
test22.uif 
test2S .uif 
test24 .uif 
test2S .uif 
test26.uif 
testSl .uif 
testSS.uif 
testS4.uif 


coded  exchange  test 

current  and  d  array 

transmission  line  test/constant  d 

transmission  line  test /travelling  wave 

3-D  coded  current  array 

3-D  field  test 

particle  mover  test 

particle  mover  test 

periodic  mover  test 

cyclotron  orbit 

cyclotron  with  E  X  B 

current  assignment  check 

short  MITL  test 

longer  MITL 

Further  E  x  B  test 

1  cavity  device/no  particle 

5  cavity  MILO/few  particles 

4  cavity  MILO/50  steps 

Symmetry  be  EM  transmission  line 

test  error  in  G2LMAT 

particle  reflection  be  test 

uniform  d  start  5  cavity  milo 

1  block  MILO  test 

20  block  MILO  test/.5000  steps 

20  block/finer  mesh  MILO  15,000  steps 

04  block  MILO 

3-D  version  of  test  1 

3-D  version  of  test  3 

3-D  version  of  test  4 


Input  dataset  testnn.uif  has  a  corresponding  output  file  o_testnnpl, 
GHOST  graphical  output  file  g.testnnpl  and  restart  file  r_testnnpl,  and 
a  PIC3D  input  file  testnn.dat.  (The  suffix  ‘pi’  is  for  a  one  processor  run. 


Task  1  Final  Report 


17 


Output  from  m  processor  runs  on  the  iPSC  have  suffices  pmtr  for  the  rth  pro¬ 
cessor  output  from  an  m  processor  hypercube).  The  runs  exercise  the  various 
libraries  of  program  units,  and  the  testnn.dat  files  provide  tests  of  the  data 
interfacing  and  of  the  main  simulation  code  PIC3D. 

Listings  of  these  test  runs  are  not  included  in  this  document  because  of  their 
length.  The  input  datasets  are  provided,  and  output  can  be  compared  to  test 
datasets  provided  with  the  report  [4]. 

In  order  to  allow  the  development  of  the  body  fitted  software  to  proceed  in 
parallel  with  the  implementation  and  testing  on  the  parallel  machines,  a  simple 
scalable  test  problem  which  had  the  same  message  passing  properties  as  the  full 
code,  but  did  not  require  the  general  geometry  metrics  and  boundary  conditions 
was  devised.  The  representative  problem  chosen  was  a  multiblock  implementa¬ 
tion  of  a  triply  periodic  electron  plasma  (see  Annex  A).  The  prototypical  test 
case  for  this,  test40.uif,  is  shown  in  Appendix  J,  together  with  the  graphical 
output  g_test40pl.  The  scatter  plot  shows  the  superposition,  viewed  along  a 
principal  mesh  axis,  of  the  superparticle  electron  positions  for  five  steps.  The 
plot  shows  the  thermal  spread  of  particles  moving  from  their  initial  positions 
on  a  lattice  of  2  x  2  x  2  particles  per  element.  For  this  particular  instance, 
the  element  net  comprised  of  a  2  x  2  x  2  cube  of  blocks,  each  block  containing 
4x4x4  elements. 

5.8  PIC3D 

The  initial  testing  of  the  data  interfacing,  control  structure  and  electromagnetic 
modules  of  PIC.3D  uses  the  .dat  files  generated  by  running  MIAIDPIC.  A  list  of 
these  tests  is  given  in  the  previous  Section.  In  addition,  further  tests  using  non- 
cartesian  element  nets  have  been  performed;  some  success  has  been  obtained 
with  cylindrical  nets  and  general  element  nets.  Cylindrical  geometry  test  which 
have  produced  reasonable  results  to  date  are  the  following: 

Transmission  Line  Tests  A  transmission  line  test  is  contained  in  the  file 
tl000.dat.  The  coaxial  transmission  line  is  2  mm  long  and  has  an  inner  radius 
of  0.5  mm  and  an  outer  radius  of  1  mm  and  the  calculation  is  run  for  0.015 
ns.  A  5x24x20  finite  element  net  is  used  and  the  calculation  is  run  for  297 
steps.  There  is  an  applied  field  boundary  at  the  input  end  and  a  resistive  wall 
boundary  with  impedance  Zq  at  the  output  end.  A  voltage  is  applied  to  the 
input  end  which  is  ramped  up  from  zero  over  0.0015  ns  and  then  held  constant. 
Figure  .3  shows  the  voltage  across  the  line  versus  time  at  the  input  and  at  the 
centre  of  the  line.  It  shows  the  correct  behaviour,  namely  that  the  signal  travels 
along  the  transmission  line  at  the  speed  of  light  and  that  there  is  no  reflection 
at  the  far  end.  The  overshoot  and  oscillation  on  the  signal  is  due  dispersion  in 
the  numerical  scheme  because  of  the  rapid  rise  time  of  the  input  signal  and  the 
coarseness  of  the  finite  element  net. 

Resonant  Cavity  Tests  A  series  of  resonant  cavity  tests  has  been  performed 
for  a  cavity  which  is  1  mm  in  radius  and  2  mm  long.  The  central  region  of  the 
device  consists  of  two  4x12x4  blocks  which  are  joined  end  to  end  each  of  which  is 
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Figure  3:  Transmission  Line  Test 


surrounded  by  three  4x4x4  blocks.  There  are  a  total  of  8  blocks  and  52  patches 
for  block  connectivity  and  boundary  conditions. 

In  these  calculations  the  magnetic  field  is  initialised  to  distributions  for 
TElOl,  TElll,  TMlOl  and  TMlll  modes  and  the  subsequent  evolution  of 
the  fields  examined  to  see  if  the  code  is  behaving  correctly.  The  initial  field  is 
computed  by  evaluating  the  analytic  form  of  the  vector  potential  at  the  centres 
of  edges  of  elements.  The  fields  should  oscillate  at  the  frequency  predicted  by 
analytic  theory.  The  tests  test  that  the  field  solver  is  using  the  correct  metrics 
on  the  axis  and  that  the  azimuthal  filtering  near  the  axis  is  working  correctly. 
The  datasets  for  these  calculations  tel01.dat,  telll.dat,  tmlOl.dat.  and 
tmlll.dat  are  in  the  TEST  subdirectory  below  P1C3D. 

The  axisymmetric  cases  TElOl  and  TMlOl  oscillate  at  a  single  frequencies 
of  196.6  GHz  and  139.4  GHz  respectively.  The  predicted  frequencies  are  197.7 
GHz  and  137.2  GHz  respectively.  The  errors  in  the  frequencies  are  approxi¬ 
mately  0.6  %  and  1.6  %  respectively. 

The  non-axisymmetric  cases  TElll  and  Thllll  oscillate  at  a  single  fre¬ 
quency  except  for  the  elements  near  the  axis.  The  measured  frequencies  are 
114.9  GHz  and  199.8  GHz  compared  with  the  theoretical  values  of  115.6  GHz 
and  197.7  GHz.  The  errors  are  0.6  %  and  1.1  %  respectively.  In  the  elements 
next  to  the  axis  the  TElll  mode  has  a  spurious  oscillation  at  about  880  GHz 
which  has  an  amplitude  of  about  50  %  of  the  signal  at  196.6  GHz  but  is  only 
significant  near  the  axis.  The  TMlll  mode  has  a  spurious  oscillation  at  about 
860  GHz  which  has  an  amplitude  of  about  50  %  of  the  signal  at  139.4  GHz  but 
is  only  significant  near  the  axis. 

The  finite  element  net  has  been  refined  so  that  there  are  24  elements  in 
the  azimuthal  direction  in  the  central  blocks  and  8  elements  in  the  azimuthal 
direction  in  the  outer  blocks.  The  numbers  of  elements  in  the  radial  and  axial 
directions  are  the  same  as  before.  The  TElll  case  gives  much  better  results 
for  the  refined  finite  element  net.  The  frequency  of  the  main  oscillation  for  the 
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TElll  calculation  is  115.5  GHz  (error  0.04  %)  and  the  spurious  oscillation  near 
the  axis  has  an  amplitude  of  about  7  %  of  the  main  oscillation.  The  TMlll 
calculations  on  the  refined  finite  element  net  are  not  significantly  different  from 
the  calculations  on  the  coarse  finite  element  net. 

When  the  finite  element  net  is  refined  by  a  factor  of  two  in  the  radial  and 
axial  direction  the  TMlll  case  oscillates  at  198.5  GHz  (error  0.4  %).  The 
spurious  oscillation  near  the  axis  is  at  1700  GHz  and  with  an  amplitude  which 
is  85  %  of  the  main  oscillation.  The  discrepancies  in  the  solutions  near  to  the 
axis  are  larger  than  expected  and  may  be  due  to  some  coding  errors  or  to  a 
non-optimal  choice  of  metrics  on  the  axis.  Further  work  is  required  to  resolve 
this  issue. 


Magnetostatic  Field  Tests  The  magnetostatic  field  tests  are  similar  to  the 
resonant  cavity  tests  except  that  the  initial  field  corresponds  to  a  magnetostatic 
solution  of  Maxwell’s  equations.  The  case  with  a  uniform  magnetic  field  in 
the  axial  direction  gives  a  field  which  remains  constant  and  uniform  to  within 
rounding  error.  The  cases  with  a  uniform  field  in  a  transverse  direction  oscillate 
around  the  steady  state  value. 


Uniform  ntacnetotiotlc  field  In  t  direction 
rediel  mefoetic  Inteoeity 


Km  ttWiU  fywwM*  MeeM 


Figure  4:  Magnetostatic  Field  Test.  The  plot  was  made  from  the  output  .tsd 
file  using  MPICTIM 

Figure  4  shows  the  evolution  of  a  the  radial  component  of  the  magnetic 
field  at  three  different  radial  locations  in  a  calculation.  The  calculation  uses  a 
8x48x8  finite  element  net  and  the  radial  field  is  plotted  at  the  centre  of  elements 
next  to  the  axis,  two  elements  away  from  the  axis  and  six  elements  away  from 
the  axis.  The  discrepancies  in  the  field  are  greatest  near  that  axis  and  have 
dropped  to  a  negligible  level  six  elements  away  from  the  axis. 


General  Geometry  Tests  General  geometry  tests  include  the  case  of  a  uni¬ 
form  magnetic  field  aligned  along  a  rectangular  cavity,  meshed  as  the  example 
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Figure  5:  The  rectangular  cavity  meshing  for  tlie  general  geometry  test 


uni20  of  Section  1.  Figure  5  illustrates  the  meshing  used.  The  field  should  re¬ 
main  invariant,  and  does  so  to  within  a  relative  error  of  5  x  10“'*.  Appendix  K 
lists  the  components  d'  and  6'  after  50  timesteps.  Fields  aligned  in  each  of  the 
other  two  coordinate  directions  have  also  been  verified  to  remain  approximately 
unchanged. 
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C  IDATA(2)  INPUT  INTEGER  DATA  ON  CHANNEL  NREAD  U.46 

C  LDATA(2)  INPUT  LOGICAL  DATA  ON  CHANNEL  NREAD  U.47 
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C  ALMN  get  A  for  Imn  mode  1.29 

C  MUNIFM  compute  A  for  uniform  magnetic  field  1.30 

C 

CL  CALCULATION  CLASS  2 

C  AMPERE  Compute  displacement  current  2.20 

C  FARADA  Advance  magnetic  field  one  step  2.21 

C  GBTOH  Compute  H  =  Gb  2.22 

C  GDTOE  Compute  E  =  Gd  2.23 

C  BCOPAT  apply  orthogonal  be  on  d  on  patches  2.25 

C  ADDVEC  add  two  vectors  in  block  together  2.26 

C  RETARD  Retard  B  by  half  a  step  2.27 

C  BCONE  apply  orthogonal  be  on  E  on  patches  2.28 

C  GLUEIO  transfer  data  to/from  gluepatch  buffers  2.30 

C  CPYVEC  copy  block  vector  2.34 

C  AVEVEC  average  two  vectors  in  block  2.35 

C  BCAXI  condition  on  d  and  E  on  axis  2.36 

C  VECFFT  perform  FFT  on  field  structures  2.50 

C  SET77  FFT  routine  2.51 

C  FFT77  FFT  routine  2.52 

C  QPASSM  FFT  routine  2.53 

C  RPASSM  FFT  routine  2.54 

C  FILTR  filter  spectral  components  2.55 

C  FILTRM  filter  spectral  components  {3D)  2.56 

C  REPAD  reset  periodic  padding  values  2.57 

C  REPADM  reset  periodic  padding  values  (3D)  2.58 

C  AMPINV  ampere  equation  in  vacuum  2.60 

C  AMPINM  ampere  equation  in  medium,  2.61 

C  AMPMAV  ampere  equation  in  vacumm  (3D)  2.62 

C  AMPMAM  ampere  equation  in  medium  (3D)  2.63 

C  ADD3D  add  D  vector  fields  {3D)  2.64 

C  AMP2DV  ampere  equation  in  vacumm  {2D)  2.65 

AMP2DM  ampere  equation  in  medium  {2D)  2.66 

OUTPUT  CLASS  3 

OP3VEC  print  vector  3.5 

OPBUF  particle  patch  exchange  buffer  op  3.14 
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C - 

CL 

c 

C  VERSION  1.00.00 

c 

CL 

C  GETMET 
C  NOPHY 
C 

CL 

C  BCSURD 
C  BCSURE 
C  BCSYM 
C  BLKIO 
C 

CL 

C  LINEGR 
C  LINEG 
C 
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PROLOGUE  CLASS  1 

compute  block  metrics 

return  block  index  corresponding  to  BP  name 

CALCULATION  CLASS  2 

apply  be  on  d  at  surfaces 
apply  be  on  E  at  surfaces 
apply  symmetry  be 

copy  data  between  block  and  gluepatch  buffers 


OUTPUT 

plot  line  graph  of  net  field 
plot  curve  for  block 


CLASS  3 


1.20 

1.21 


2.24 

2.27 

2.29 

2.32 


3.13 

3.14 


o  o 
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c - 

CL 

C 

C  VERSION 
C 

CL 

C  XACMAT 
C 

CL 

C  SETCUR 
C  MOVCUR 
C  ACCEL 
C  EMITEL 
C  BEAMIN 
C  PPSORT 
C  PCXFRM 
C  ROTSPV 
C  QSHARE 
C  ASSCUR 
C  MERGE 
ACCEL2 
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PROLOGUE  CLASS  1 

find  acceleration  transformation  matrices 

CALCULATION  CLASS  2 

Initialise  current  arrays 
Move  particles  and  compute  currents 
Update  particle  m.omenta 
Emit  electrons  from  cathode  surfaces 
Inject  electron  beam 
sort  particles  into  patch  buffer 
transform  coord  to  target  block  coords 
rotate  and  scale  p  and  v  to  new  x 
assign  charge  to  block  mesh 
assign  current  from  trajectory  segment 
merge  2  ascending  lists  of  coordinates 
Update  particle  momenta  for  cart/polar  only 


1.15 


2.10 

2.11 

2.12 

2.13 

2.14 

2.15 

2.16 

2.17 

2.18 

2.19 

2.20 
2.21 
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c - 

CL 

c 

C  VERSION  1.00.00 

c 

CL 

C  PINIT 
C  SETPAR 
C 

CL 

C  PART 10 
C  PBLINK 
C  INJECT 
C 

CL 

C  PARTOP 
C  SCAPLT 
C  PARSUM 
C 
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PROLOGUE  CLASS  1 

initialise  particles 
set  particle  addressing 

CALCULATION  CLASS  2 

load  particles  into/from  exchange  buffers 
manage  particle  buffering  between  blocks 
apply  particle  injection  be 

OUTPUT  CLASS  3 

particle  table  summary 
particle  scatter  plot 
summary  of  particle  storage 


1.13 

1.14 


2.14 

2.17 

2.19 


3.6 

3.7 
3.9 
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c - 

CL  INDEX  OF  SUBPROGRAMS 

C 

C  VERSION  1.00.00  JWE/NJB  Culham  Laboratory  Feb  199« 
C 


CL 

CONTROL 

CLASS 

0 

C 

p 

RECORD 

read  and  write  restart  data 

0.5 

CL 

PROLOGUE 

CLASS 

1 

C 

CLEAR 

set  common  variables  and  arrays 

1.2 

c 

PRESET 

Set  default  values 

1.3 

c 

RESUME 

resume  from  previous  run 

1.7 

c 

MOUDAT 

output  IV  from  MIMDPIC  for  PIC3D 

1.9 

c 

ADRDIM 

set  addressing  and  dimensioning  parameters 

1.39 

c 

SETADP 

set  addressing  parameters 

1.44 

c 

SETDIP 

set  dimensioning  parameters 

1.45 

c 

SETFDP 

set  fixed  dimensioning  parameters 

1.46 

c 

SETPCN 

set  physical  constants 

1.47 

c 

Q 

DATA 

Input  data  specific  to  run 

1.4R 

CL 

OUTPUT 

CLASS 

3 

c 

OUTDAT 

write  out  data  file 

3.1 

c 

FRAMED 

label  and  advance  frame 

3.10 

c 

GHINIT 

initialise  graphical  output 

3.11 

C 
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VERSION  1.00.00 


C 

CL 

C 

C 

C 

CL 

C 

C 

C 

CL 

C 

c 

CL 

c 

c 


SETDMP 

NODASG 


XPATCH 

XPART 


PATOP 


CSEND 
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PROLOGUE  CLASS  1 

setup  tables  for  xpatch 
assign  processes  to  processors (nodes) 

CALCULATION  CLASS  2 

exchange  gluepatch  buffers 
exchange  particle  coordinates 

OUTPUT  CLASS  3 

patch  diagnostic  output 

UTILITIES  CLASS  P 

dummies  for  Intel  routines 


1.11 

1.12 


2.31 

2.33 
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CL 

INDEX  OF  SUBPROGRAMS 

C  VERSION  1 

N  J  Brealey  AEA  Technology 

June 

1994 

CL 

CALCULATION 

CLASS 

2 

C 

DIAINT 

initialise  diagnostics 

2.18 

c 

DIACOL 

diagnostic  collection 

2.19 

c 

DIAOUT 

output  diagnostic  data 

2.20 

c 

DIAFIN 

finish  diagnostic  output 

2.21 

CL 

OUTPUT 

CLASS 

3 

c 

TSDHED 

read  or  write  TSD  header 

3.1 

c 

TSDENT 

read  or  write  TSD  entry 

3.2 

c 

TSDSET 

setup  TSD  output 

3.3 

c 

TSFLAG 

set  flags  for  time  set 

3.4 

c 

GFILNM 

get  file  name 

3.5 

c 

SETDWS 

set  diagnostic  workspace  origins 

3.6 

c 

GRLSET 

setup  GRL  output 

3.7 

c 

TSDOUT 

output  TSD  data 

3.8 

c 

GRLOUT 

output  GRL  data 

3.9 

c 

GRLPLT 

plot  GRL  data 

3.10 

c 

LPLOT 

plot  a  set  of  lines 

3.11 

CL 

EPILOGUE 

CLASS 

4 

C 

lEDL 

line  integral  of  E 

4.1 

c 

lEDLM 

line  integral  of  E  (3D) 

4.2 

c 

IHDL 

line  integral  of  H 

4.3 

c 

IHDLM 

line  integral  of  H  (30) 

4.4 

c 

IDDS 

surface  integral  of  D 

4.5 

c 

IDDSM 

line  integral  of  E  (3D) 

4  .  6 

c 

IBDS 

surface  integral  of  D 

4.7 

c 

IBDSM 

surface  integral  of  D  (3D) 

CO 

c 

IBHV 

volume  integral  of  B.H 

4.9 

c 

IBHVM 

volume  integral  of  B.H  (3D) 

4.10 

c 

lEDV 

volume  integral  of  E.D 

4.11 

c 

lEDVM 

volume  integral  of  E.D  (3D) 

4.12 

c 

ARCLN 

compute  arc  lengths 

4.13 

c 

ARCLNM 

compute  arc  lengths  (3D) 

4.14 

c 

E123 

physical  component  of  E 

4.15 

c 

E123M 

physical  component  of  E  (3D) 

4.16 

c 

H123 

physical  component  of  H 

4.17 

c 

H123M 

physical  component  of  H  (3D) 

4.18 

c 

SUCLR 

clear  SU  workspace 

4.19 

c 

SUSCL 

scale  SU  workspace 

4.20 

c 

SUADD 

add  SU  workspace 

4.21 

c 

SUMAX 

max  and  min  of  SU  workspace 

4.22 

c 
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c - 

CL  INDEX  OF  SUBPROGRAMS 

C 

C  VERSION  1.00.00  WA  Culham  Laboratory  Mar  1994 
C 

CL  CONTROL  CLASS  0 

RECORD  read  and  write  restart  data  0.5 

PROLOGUE  CLASS  1 

LABRUN  Label  the  run  1.1 

CLEAR  set  common  variables  and  arrays  1.2 

PRESET  Set  default  values  1.3 

AUXVAL  SET  AUXILLIARY  VALUES  1.5 

START  read  records  1.8 

RDFMTl  read  tagged  format  1  1.11 

RDFMT2  read  tagged  format  2  1.12 

RDFMT3  read  tagged  format  3  1.13 

RDFMT4  read  tagged  format  A  1.14 

RDFMT5  read  tagged  format  5  1.15 

RDFMT6  read  tagged  format  6  1.16 

RDFMT7  read  tagged  format  1  1.17 

PREORD  cube  geometry  ordering  1.31 

PRECLA  set  up  class  list  1.32 

PREQUD  quadrature  for  metric  1.33 

GETECW  compute  block  basis  vectors  1.39 

SCALEV  scale  covariant  vectors  in  block  1.40 

DATA  Input  data  specific  to  run  1.4R 

CL  CALCULATION  CLASS  2 

C  STEPON  control  computation  2.1 

C  VCLAS  validate  classes  2.3 

C  SEXION  sectional  information  2.4 

C  INCODE  input  conversion  and  sorting  2.5 

C  VCOVER  check  all  blocks  covered  2.6 

C  VORI  check  relative  orientations  2.7 

C  VCONN  check  block  connectivity  2.8 

C  CPATBK  surface  patch  -  block  tables  2.12 

C  CLIPAT  line  patches  2.13 

C  PATEDG  sort  patches  by  edge  2.14 

C  PTRAIL  patch  trails  2.15 

C  PEGGY  peg  into  main  line  notation  2.16 

C  CGLOBV  global  positioning  variables  2.17 

C  CGPHYA  global  physical  attribute  table  2.18 

C  GETSP  compute  boundary  condition  attributes  2.19 

C  NAMLIS  namelist  handling  2.20 

C  NAMINP  real  namelist  input  2.21 

C  OXCONV  CONVERT  OX  NOTATION  TO  INTERNAL  REPN  2.22 

C  LOCPOS  local  positions  2.23 

C  XYZROT  block  vectors  and  rotation  matrices  2.24 

C  CROTBK  ultimate  global  system  2.25 

C  P2LMAT  provisional  global  to  local  coord  rotation  mat  2.26 

C  GETCDT  calculate  timestep  *  (speed  of  light)  2.27 

C  CALCDT  calculate  timestep  *  c  for  block  2.28 

C  GETSCA  scaling  to  dimensionless  variables  2.29 

C 

CL  OUTPUT  CLASS  3 

C  OUTPUT  control  output  3.1 

C  DIAGN  Diagnostic  output  3.3 

C 

CL  EPILOGUE  CLASS  4 

C  TESEND  test  for  completion  of  run  4.1 

C  ENDRUN  terminate  the  run  4.2 
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CL  INDEX  OF  SUBPROGRAMS 

C 

C  VERSION  1.00.00  N  J  Brealey  AEA  Technology  Feb  199A 

C 


CL 

CONTROL 

CLASS 

0 

C 

r- 

COTROL 

Control  the  run 

0.3 

L. 

CL 

PROLOGUE 

CLASS 

1 

c 

LABRUN 

label  the  run 

1.1 

c 

AUXVAL 

set  auxilary  values 

1.5 

c 

INITAL 

setup  physical  initial  conditions 

1.6 

c 

RESTOR 

restore  device  description 

1.9 

c 

STRAGE 

allocate  storage 

1.10 

c 

MIMDAT 

read  data  for  MIMDPIC  commons 

1.11 

c 

BINIT 

compute  initial  B 

1.12 

c 

r' 

MINDAT 

I.V.  input  data  from.  MIKDPIC 

1.13 

L. 

CL 

CALCULATION 

CLASS 

2 

C 

STEPON 

advance  calculation  one  step 

2.1 

c 

MOVEl 

inject  and  particle  move  1 

2.2 

c 

DFILTR 

filter  Ddot 

2.3 

c 

MOVE2 

particle  move  2 

2. A 

c 

CFILTR 

filter  current 

2.5 

c 

CSPLIT 

split  current  for  subcycling 

2.6 

c 

DDOT 

compute  ddot 

2.1 

c 

BUFUNL 

unload  glue  patch  buffers 

2.8 

c 

INJECT 

inject  particles 

2.9 

c 

DTOE 

convert  D  to  E  in  blocks 

2.10 

c 

BTOH 

convert  B  to  K 

2.11 

c 

ADVB 

advance  B 

2.12 

c 

RETB 

retard  B  half  a  step 

2.13 

c 

NEWP 

update  mom.cnta 

2.1A 

c 

EBEXT 

add  external  field 

2.15 

c 

r' 

DDOT2 

alternative  version  od  DDOT/BTOH 

2.18 

CL 

OUTPUT 

CLASS 

3 

C 

OUTPUT 

control  the  output 

3.1 

C 

MIMOUT 

output  for  KIKDPIC 

3.2 

C 

RESTRT 

write  to  disc 

3.-) 

L. 

CL 

EPILOGUE 

CLASS 

A 

C 

TESEND 

test  for  completion  of  run 

k.2 

C 
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CL  INDEX  OF  SUBPROGRAMS 

C 

C  VERSION  1.00.00  JWE  Culham  Laboratory  Feb  1994 
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CL 

CONTROL 

CLASS  0 

C 

BASIC 

Version  without  input  file 

read  for  MIMD 

0.1 

c 

MODIFY 

Version  to  create  separate 

ouput  files 

0.2 

c 

COTROL 

Control  the  run 

0.3 

c 

EXPERT 

test  data  dump  version 

0.4 

c 

Q 

RECORD 

Read  and  write  restart  coordinates 

0.5 

CL 

PROLOGUE 

CLASS  1 

c 

LABRUN 

Label  the  run 

1.1 

c 

CLEAR 

Clear  common  arrays  and  variables 

1.2 

c 

PRESET 

Set  default  values,  not  MYPRES  etc 

1.3 

c 

DATOLD 

old  version  of  data  cls4r 

1.4 

c 

AUXVAL 

Compute  auxiliary  values 

1.5 

c 

INITAL 

Define  geometrical  and  physical  IV 

1.6 

c 

RESUME 

resume  from  previous  run 

1.7 

c 

START 

start  the  calculation 

1.8 

c 

Q 

MOUDAT 

output  IV  for  PIC3D 

1.9 

CL 

CALCULATION 

CLASS  2 

c 

Q 

STEPON 

step  on  the  calculation 

2.1 

CL 

OUTPUT 

CLASS  3 

C 

f' 

OUTPUT 

Control  the  output 

3.1 

CL 

epilogue" 

CLASS  4 

C 

TESEND 

Test  for  completion  of  run 

4.1 

c 

r* 

ENDRUN 

Terminate  the  run 

4.2 

CL 

DIAGNOSTICS 

CLASS  5 

C 

CLIST 

Dump  common  variables 

5.2 

C 

ARRAYS 

Dump  common  arrays 

5.3 

CL 

UTILITIES 

CLASS  X 

C 

GLBCLR 

clear  global  array  quantities 

X.l 

c 

RECTBK 

set  up  rectangular  brick  blocktype 

X.l 

c 

GLBCLI 

Dump  common  variables 

X.2 

c 

GRDPLT 

plot  parallel  projection  of  block  grid 

X.2 

c 

GLBARR 

Dump  common  arrays 

X.3 

c 

NETPLT 

plot  finite  element  net 

X.3 

c 
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B  Sample  preprocessor  input  dataset 


pbc.uif 

o_pbcl 

'NLRES  L  False  for  NEWRUN,True  for  RESET'  F 
pbcl 

Periodic  section  of  rectangular  waveguide 
using  two  cartesian  blocks 
conducting  walls  in  1  and  3  directions 
periodic  be  in  2  direction 


'CHLAB5 

A 

**Label  available  to  programmer 

1.1 

/ 

/ 

'NDIARY 

I 

**Channel 

for  diary 

1.2 

r 

/ 

'NIN 

I 

**Current 

input  channel 

1.2 

$ 

/ 

'NLEDGE 

I 

**Channel 

for  restart  records 

1.2 

t 

/ 

'NONLIN 

I 

**Channel 

for  input-output 

1.2 

t 

/ 

'NOUT 

I 

**Current 

output  channel 

1.2 

t 

/ 

'NPRINT 

I 

**Channel 

for  printed  output 

1.2 

/ 

/ 

'NPUNCH 

I 

**Channel 

for  card  output  (or  equivalent) 

1.2 

/ 

/ 

'NREC 

I 

**Current 

record  number 

1.2 

/ 

/ 

'NRUN 

I 

**Maximum 

number  of  steps 

1.2 

t 

/ 

' NADUMP 

lA 

**Codes  for  array  dumps 

1.9 

r 

/ 

'NPDUMP 

lA 

**Codes  for  dumping  points 

1.9 

t 

/ 

' NVDUMP 

lA 

**Codes  for  dumping  arrays 

1.9 

t 

/ 

'NLCHED 

L 

**.TRUE. 

if  class  0  report-head  required 

1.9 

/ 

/ 

'NLHEAD 

LA 

**.TRUE. 

if  class  1-9  report-heads  required 

1.9 

t 

/ 

'NLOMTl 

LA 

**Class  1 

subprogram  selector 

1.9 

t 

/ 

'NLOMT2 

LA 

**Class  2 

subprogram  selector 

1.9 

t 

/ 

'NLOMT3 

LA 

**Class  3 

subprogram  selector 

1.9 

/ 

/ 

'NLREPT 

L 

**.TRUE. 

if  report  required 

1.9 

/ 

/ 

BL  blockl 

cube  uniforml 

BL  block2 

cube  uniforml 

BG  cube  regular  cubic 

lattice  -6  -10  -16 

SCUBREG  RXMAX=1 . ,  RyMAX=2 . , RZMAX=3 . / 
BP  uniforml  uniform 
SUNIFRM  EPSR=1.,RMUR=1./ 

PP  pcond  pcond 

PA  pcond  sameas  pcond 

BC  blockl 

W=pcond 

S=blockl  :N 

N=blockl:S 

E=block2 :W 

U=pcond 

Down  =poond 

EN 

BC  block2 
W=blockl :E 
S=block2 :N 
N=block2 :S 
E=pcond 
U=pcond 
D=pcond 
EN 

OR  blockl  block2 


'  94/09/22 
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C  PIC3D  input  dataset  template 


GLBLIB/.NEWRUN_DAT 


94,'m'I6 

16:09:06 


o  newrun 

'NLRES  L 

expOOOOl 

CHLABl 

CHLAB2 

CHLAB3 

CHLAB4 

' CHLAB5 

False  for  NEWRUN, True  for  RESET'  / 

8  character  label  for  run  sequence 
replace  these  four  lines 
by  run  labelling  information 
which  will  appear  at  the  start 
of  the  NOUT  channel  output 

A  **Label  available  to  programmer 

1.1 

/ 

/ 

'NDIARY 

I 

**Channel  for  diary 

1.2 

/ 

/ 

'NIN 

I 

**Current  input  channel 

1.2 

t 

/ 

' NLEDGE 

I 

**Channel  for  restart  records 

1.2 

/ 

/ 

'NONLIN 

I 

**Channel  for  input-output 

1.2 

t 

/ 

'NOUT 

I 

**Current  output  channel 

1.2 

0 

/ 

'NPRINT 

I 

**Channel  for  printed  output 

1.2 

0 

/ 

'NPUNCH 

I 

**Channel  for  card  output  (or  equivalent) 

1.2 

0 

/ 

'NREC 

I 

**Current  record  number 

1.2 

0 

/ 

'NRUN 

I 

**Maximum  number  of  steps 

1.2 

0 

/ 

' NADUMP 

lA 

**Codes  for  array  dumps 

1.9 

0 

/ 

'NPDUMP 

lA 

**Codes  for  dumping  points 

1.9 

0 

/ 

' NVDUMP 

lA 

**Codes  for  dumping  arrays 

1.9 

0 

/ 

'NLCHED 

L 

**.TRUE.  if  class  0  report-head  required 

1.9 

0 

/ 

'NLHEAD 

LA 

**.TRUE.  if  class  1-9  report-heads  required 

1.9 

0 

/ 

'NLOMTl 

LA 

**Class  1  subprogram  selector 

1.9 

0 

/ 

'NLOMT2 

LA 

**Class  2  subprogram  selector 

1.9 

0 

/ 

'NLOMT3 

LA 

**Class  3  subprogram  selector 

1.9 

t 

/ 

'NLREPT 

L 

**.TRUE.  if  report  required 

1.9 

0 

/ 

'COUR 

R 

*max  courant  number 

3.1 

0 

/ 

'DT 

R 

*Timestep 

3.1 

0 

/ 

' SCALE 

RA 

*scale  to  SI  factors 

3.1 

0 

/ 

'NSUBS 

I 

*number  of  em  steps  per  particle  step 

3.1 

0 

/ 

'ROTG2L 

RA 

"block  global  to  local  coord  rotation 

4  ;i 

0 

/ 

' XYZBLK 

RA 

"block  corner  global  coords 

4.1 

0 

/ 

' MBKPAT 

lA 

"block  to  patch  hoc  table 

4.1 

0 

/ 

' MBKTYP 

lA 

"block  to  blocktype  pointer 

4.1 

0 

/ 

' MEMBCA 

lA 

"EM  patch  to  BC  attribute  table 

4.1 

0 

/ 

' MPABCA 

lA 

"particle  patch  to  BC  attribute  table 

4.1 

0 

/ 

'MPATBK 

lA 

"patch  to  block  pointer 

4.1 

0 

/ 

' MPATO 

lA 

"location  code  of  patch  origin  on  blocks 

4.1 

0 

/ 

' MPATX 

lA 

"location  code  of  patch  extreme  on  blocks 

4.1 

0 

/ 

'MPATYP 

lA 

"patch  type  index 

4.1 

0 

/ 

' NBLOCK 

I 

"number  of  blocks 

4.1 

0 

/ 

'NPATCH 

I 

"number  of  patches 

4.1 

0 

/ 

'NXTPAT 

lA 

"link  to  next  patch  on  block 

4.1 

0 

/ 

' LBLAS 

lA 

"list  of  block  addressing  structures 

4.2 

0 

/ 

' LECOVA 

lA 

"list  of  ECOV  Addressing 

4.2 

0 

/ 

' LGBHAD 

lA 

"list  of  G  addressing  for  H=G  b 

4.2 

0 

/ 

' LGDEAD 

lA 

"list  of  G  addressing  for  E=G  d 

4.2 

0 

/ 

' LOAGBH 

lA 

"location  of  GBH  addr.  in  LGBHAD 

4.2 

0 

/ 

' LOAGDE 

lA 

"location  of  GDE  addr.  in  LGDEAD 

4.2 

0 

/ 

'  LOBLAS 

lA 

"location  of  origins  of  block  addresses 

4.2 

0 

/ 

' LOECOA 

lA 

"location  of  addressing  in  LECOVA 

4.2 

0 

/ 

' LOECOV 

lA 

"location  of  origins  in  ECOV  array 

4.2 

0 

/ 

' LORGBH 

lA 

"location  of  origins  of  GBH  in  G  array 

4.2 

0 

/ 

'  LORGDE 

lA 

"location  of  origins  of  GDE  in  G  array 

4.2 

0 

/ 

'  MBPROP 

lA 

"block  properties 

4.2 

0 

/ 

'NBTYPE 

I 

"number  of  block  types 

4.2 

0 

/ 

' SPATR 

RA 

"array  of  particle  species  attributes 

4.3 

0 

/ 

' LOCOOR 

lA 

""Location  of  Origins  of  COORdinates 

4.3 

0 

/ 

' LOPARA 

lA 

"Location  of  Origin  of  blocks  in  LPARAS 

4.3 

0 

/ 

' LPARAS 

lA 

"List  of  addressing  for  coordinates 

4.3 

0 

/ 

'NSPEC 

I 

"number  of  particle  species 

4.3 

0 

/ 

' MBKPES 

lA 

"block  to  processor  pointer 

4.4 

0 

/ 

' MPRBLK 

lA 

"process  to  block  hoc  table 

4.4 

0 

/ 

'MSTAT 

lA 

"Status  of  message  array 

4.4 

0 

/ 

'  NXTBLK 

lA 

"next  block  in  processor 

4.4 

0 

/ 

'  MDBCA 

I 

"Dimension  of  be  attribute  arrays 

4.6 

0 

/ 

'MDIOB 

I 

""Dimension  of  10  Buffer  arrays 

4.6 

0 

/ 

' MDPART 

I 

""Dimension  of  PARTicle  coordinate  arrays 

4.6 

0 

/ 

'MECDIM 

I 

"ECov  array  DIMension 

4.6 

0 

/ 

'MGDIM 

I 

"Metric  G  array  DIMension 

4.6 

0 

/ 

'NDSPS 

lA 

""numbers  of  domain  sets  in  each  plot  set 

5.0 

0 

/ 

'NFSPS 

lA 

""numbers  of  field  sets  in  each  plot  set 

5.0 

0 

/ 

'NPS 

I 

""number  of  plot  sets 

5.0 

0 

/ 

'NTSPS 

lA 

""numbers  of  time  sets  in  each  plot  set 

5.0 

0 

1 

'CFMTPS 

AA 

""file  format  for  each  plot  set 

5.1 

0 

/ 

'MFLFS 

lA 

""number  of  fields  in  each  field  set 

5.2 

0 

/ 

' NFLFS 

lA 

""numbers  of  fields  in  each  field  set 

5.2 

0 

/ 

'NFS 

I 

""number  of  field  sets 

5.2 

0 

/ 

94/08/16 

16:09:06 

GLBLIBANEWRUN_DAT 

'NPSFS 

lA 

**plot  sets  of  each  field  set 

5.2 

r 

/ 

'MDODS 

lA 

"‘number  of  domains  in  each  domain  set 

5.3 

t 

/ 

'NDODS 

lA 

""numbers  of  domains  in  each  domain  set 

5.3 

t 

/ 

'NDS 

I 

""number  of  domain  sets 

5.3 

/ 

/ 

'NPSDS 

lA 

""plot  sets  of  each  domain  set 

5.3 

t 

/ 

' TIMETS 

RA 

""start,  stop,  output,  average  times,  time 

z5.'1 

! 

/ 

'NPSTS 

lA 

""plot  sets  of  each  time  set 

5.4 

t 

/ 

'NSTPTS 

lA 

""start,  stop,  output,  average  steps,  step 

z5 . 4 

/ 

/ 

'NTS 

I 

""number  of  time  sets 

5.4 

/ 

/ 

' CNLTIM 

AA 

""long  names  of  time 

5.5 

t 

/ 

'CNSTIM 

AA 

""short  names  of  time 

5.5 

t 

/ 

'CNUFRQ 

AA 

""units  of  frequency 

5.5 

/ 

/ 

'CNURAT 

AA 

""units  of  rate  of  change 

5.5 

f 

/ 

'CNUTIM 

AA 

""units  of  time 

5.5 

! 

/ 

' SOLDO 

RA 

""zero  and  unit  of  arc  length 

5.6 

t 

/ 

' SLIMDO 

RA 

""domain  arc  length  limits 

5.6 

t 

/ 

'MSUDO 

lA 

""number  of  subdomiains  in  each  domain 

5.6 

t 

/ 

'NCOLDO 

lA 

""domain  colours 

5.6 

t 

/ 

'NDO 

I 

""number  of  domiains 

5.6 

t 

/ 

'NDSDO 

lA 

""domain  sets  of  each  domain 

5.6 

t 

/ 

'NSUDO 

lA 

""numbers  of  subdom.ains  in  each  domain 

5.6 

t 

/ 

'CNLDO 

AA 

""long  names  of  domains 

5.7 

t 

/ 

'CNSDO 

AA 

""short  names  of  domains 

5.7 

0 

/ 

'FLIMFL 

RA 

""field  limits 

5.8 

0 

/ 

' SCLFL 

RA 

""zero  and  unit  of  field 

5.8 

0 

/ 

'NCOLFL 

lA 

""field  colours 

5.8 

0 

/ 

'NFL 

I 

""number  of  fields 

5.8 

0 

/ 

'NFSFL 

lA 

""field  sets  of  each  field 

5.8 

0 

/ 

'CFLSPC 

AA 

""field  specification 

5.9 

0 

/ 

'CNLFL 

AA 

""long  names  of  fields 

5.9 

0 

/ 

'CNSFL 

AA 

""short  names  of  fields 

5.9 

0 

/ 

'CNUFL 

AA 

""units  of  fields 

5.9 

0 

/ 

'NPSU 

lA 

""primary  subdomain  of  subdomain 

6.0 

0 

/ 

'NSU 

I 

""number  of  subdomains 

6.0 

0 

/ 

'NSUBLK 

lA 

""subdomain  block  num.ber 

6.0 

0 

/ 

'NSUD 

lA 

""subdomain  directions 

6.0 

0 

/ 

'NSUDOM 

lA 

""subdomain  domain  number 

6.0 

0 

/ 

'NSUDS 

lA 

""subdomain  domain  set  number 

6.0 

0 

/ 

'NSUFS 

lA 

""subdomain  field  set  number 

6.0 

0 

/ 

'NSUO 

lA 

""subdomain  o  element 

6.0 

0 

/ 

'NSUPS 

lA 

""subdomain  plot  set  number 

6.0 

0 

/ 

' NSUTS 

lA 

""subdomain  time  set  number 

6.0 

0 

/ 

'NSUX 

lA 

""subdomain  x  element 

6.0 

0 

/ 

'NXTPSU 

lA 

""next  primary  subdomain 

6.0 

0 

/ 

' AMODE 

RA 

"box  size  for  modes 

6.1 

0 

/ 

'BUNI 

RA 

"uniform  initial  B 

6.1 

0 

/ 

'CMODE 

RA 

"mode  amplitudes 

6.1 

0 

/ 

'NBEXT 

I 

"use  initial  B  as  external  B 

6.1 

0 

/ 

'NINIT 

I 

"initial  field  type 

6.1 

f 

/ 

'NMODE 

lA 

"mode  numbers 

6.1 

0 

/ 
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D  CRONUS  unit  test  data 


umm 

test.dat 

o_test 

'NLRES  L  False  for  NEWRUN,True  for  RESET'  / 
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J  Eastwood,  14'^  Arter,  N  J  Brealey,  R  Hock 


E  OLYMPUS  unit  test  data 


94/09/19 

14:36:40 


os  Implementation 
Network  node: 

OS  release: 

OS  version: 
Hardware  name: 

Time:  14:36:35 
Date:  19-Sep-94 


SunOS 

rfesunh 

5.2 

Generlc_100999-51 

sun4c 


Number  of 

command  line  arguments 

Argument 

0: 

. / olytst 

Argument 

.1: 

one 

Argument 

2: 

two 

Argument 

3: 

three 

3 


Checking  CONCAT  function: 
SunOSS  .  2 


Run  time: 
User  time: 
System  time: 
Total  time: 


3.47000 

3.12000 

0.350000 

3.47000 


test.Solaris2 


48  J  W  Eastwood,  IF  Arter,  N  J  Brealey,  R  Hockney 

F  NETBLK  unit  test  data 


another  projection?  y/n 


test.dat 


wmm 

tmiM 

o_test 

'NLRES  L  False  for  NEWRUN,True  for  RESET'  / 
testOl 

NETBLK  unit  test  data 

To  run,  type  xtest  test.dat 

In  the  interactive  dialogue, 

accept  defaults  to  get  standard  test  OP. 


' CHLAB5 

A 

**Label  available  to  programmer 

1.1 

'NDIARY 

I 

**Channel 

for  diary 

1.2 

'NIN 

I 

**Current 

input  channel 

1.2 

'NLEDGE 

I 

*'*Channel 

for  restart  records 

1.2 

'NONLIN 

I 

**Channel 

for  input-output 

1.2 

'NOUT 

I 

**Current 

output  channel 

1.2 

'NPRINT 

I 

**Channel 

for  printed  output 

1.2 

'NPUNCH 

I 

**Channel 

for  card  output  (or  equivalent) 

1.2 

'NREC 

I 

‘♦Current 

record  number 

1.2 

'NRUN 

I 

“Maximum 

number  of  steps 

1.2 

'NADUMP 

lA 

“Codes  for  array  dumps 

1.9 

'NPDUMP 

lA 

“Codes  for  dumping  points 

1.9 

'NVDUMP 

lA 

“Codes  for  dumping  arrays 

1.9 

'NLCHED 

L 

“.TRUE. 

if  class  0  report-head  required 

1.9 

'NLHEAD 

LA 

“  .TRUE. 

if  class  1-9  report-heads  required 

1.9 

'NLOMTl 

LA 

“Class  1 

subprogram  selector 

1.9 

'NL0MT2 

LA 

“Class  2 

subprogram  selector 

1.9 

'NL0MT3 

LA 

“Class  3 

subprogram  selector 

1.9 

'NLREPT 

L 

“.TRUE. 

if  report  required 

1.9 

mmm 

m39m  O  test 
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Task  1  Final  Report 


H  EMBLK  unit  test  data 


i  wmm 

\  xtest.log 

rfesunh%  xtest 
??  DIR=1,  2,  or  3? 

1 

IDIR  =  1 

LBLAS  (I)  : 

24  16  8  1  25  425  3825  0  30  0 

Nl=  24  N2=  16  N3=  8 

INC1=  1  INC2=  25  INC3=  425 

IOR(l)=  1  IOR(2)=  3826  I0R(3)=  "7651 

DATAO  Values  as  set 

DATA: 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 

2.0600E+02  6.7748E+01  2.0854E+02  3.0735E401  2.1039E+02  9.3786E-01 

2.1200E+02  -2.6972E+01  2.1361E+02  -5.7078E+01  2.1546E+02  -9.4925E+01 

2.1800E+02  -1.5284E+02  2.2239E+02  -2.7473E+02  2.3439E+02  -8.7585E+02 

2.7560E+03 


ICOMP=  1 
IZ=  0 
0 

Initial  data  (first  six) : 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 

IY=  1 

Initial  data  (first  six) : 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 

IZ=  1 
IY=  0 

Initial  data  (first  six) : 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 

IY=  1 

Initial  data  (first  six): 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 

ICOMP=  2 
12=  0 
IY=  0 

Initial  data  (first  six)  : 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 

IY=  1 

Initial  data  (first  six)  : 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 

12=  1 
IY=  0 

Initial  data  (first  six)  : 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 

IY=  1 

Initial  data  (first  six) : 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 

I COMP =  3 

12=  0 
IY=  0 

Initial  data  (first  six)  : 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 

IY=  1 

Initial  data  (first  six) : 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 

12=  1 
IY=  0 

Initial  data  (first  six) : 

2.7560E+03  7.3445E+02  1.8961E+02  2.3708E+02  2.0161E+02  1.2344E+02 


xtest.log 


94/09/l‘> 

09:27:55 

IY=  1 


Initial  data  (first  six)  : 

2.7560E  +  03  7.3'!45E+02  1.8961E  +  02 

DATAO  Before  analysis 

2.3708E+02 

2.0161E+02 

1.2344E+02 

DATA: 

2.7560E+03  7.3445E+02 
2.0600E+02  6.7748E+01 
2.1200E+02  -2.6972E+01 
2.1800E+02  -1.5284E+02 
2.7560E+03 

1.8961E+02 

2.0854E+02 

2.1361E+02 

2.2239E+02 

2.3708E+02 

3.0735E+01 

-5.7078E+01 

-2.7473E+02 

2.0161E+02 

2.1039E+02 

2.1546E+02 

2.3439E+02 

1.2344E+02 

9.3786E-01 

-9.4925E+01 

-8.7585E+02 

IDIR=  1  Just  before  first  CALL  VECFFT 
Just  after  analysis 

IOR(l)=  1  IOR(2)=  3826  IOR(3)=  7651 

IDIR=  1  Just  after  first  CALL  VECFFT 

DATAO  After  analysis 

DATA: 

2.7560E+03  7.3445E+02 
2.0600E+02  6.7748E+01 
2.1200E+02  -2.6972E+01 
2.1800E+02  -1.5284E+02 
2.7560E+03 

1.8861E+02 

2.0854E+02 

2.1361E+02 

2 .2239E+02 

2.3708E+02 

3.0735E401 

-5.7078E+01 

-2.7473E+02 

2.0161E+02 

2.1039E+02 

2.1546E+02 

2.3439E+02 

1.2344E+02 

9.3786E-01 

-9.4925E+01 

-8.7585E+02 

ICOMP=  1 

IZ=  0 

IY=  0 

A(I)  after  anal  (first 
2.0000E+02  l.OlOOE+02 

IY=  1 

six)  : 

-5 .0500E+01 

1.0200E+02 

-5.1000E+01 

1.0300E+02 

Ad)  after  anal  (first 
2.0000E+02  l.OlOOE+02 

IZ=  1 

IY=  0 

six)  : 

-5.0500E+01 

1.0200E+02 

-5.1000E+01 

1.0300E+02 

A(I)  after  anal  (first 
2.0000E+02  l.OlOOE+02 

IY=  1 

six)  : 

-5.0500E+01 

1.0200E+02 

-5.1000E+01 

1.0300E+02 

A(I)  after  anal  (first 
2.0000E+02  l.OlOOE+02 

six)  : 

-5.0500E+01 

1.0200E+02 

-5.1000E+01 

1.0300E+02 

ICOMP=  2 

IZ=  0 

IY=  0 

A(I)  after  anal  (first 
2.0000E+02  l.OlOOE+02 

IY=  1 

six)  : 

-5 .0500E+01 

1.0200E+02 

-5.1000E+01 

1 .0300E  +  02 

A(I)  after  anal  (first 
2.0000E+02  l.OlOOE+02 

IZ=  1 

IY=  0 

six)  : 

-5.0500E+01 

1.0200E+02 

-5.1000E+01 

1.0300E+02 

A(I)  after  anal  (first 
2.0000E+02  l.OlOOE+02 

IY=  1 

six)  : 

-5.0500E+01 

1.0200E+02 

-5.1000E+01 

1.0300E+02 

A(I)  after  anal  (first 
2.0000E+02  l.OlOOE+02 

six)  : 

-5.0500E+01 

1.0200E+02 

-5.1000E+01 

1.0300E+02 

ICOMP=  3 

IZ=  0 

IY=  0 

Ad)  after  anal  (first 
2.0000E+02  l.OlOOE+02 

IY=  1 

six)  : 

-5.0500E+01 

1.0200E+02 

-5.1000E+01 

1 .0300E+02 

Ad)  after  anal  (first 
2.0000E+02  l.OlOOE+02 

six)  : 

-5.0500E+01 

1.0200E+02 

-5.1000E+01 

1.0300E+02 

IZ=  1 


94/09/19 
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IY= 


0 


A(I)  after  anal  (first  six): 
2.0000E+02  l.OlOOE+02  -5.0500E+01 
IY=  1 

A(I)  after  anal  (first  six): 
2.0000E+02  l.OlOOE+02  -5.0500E+01 

IDIR=  1  Just  before  second  CALL  VECFFT 
Just  after  synttiesis 

I0R{1)=  1  IOR(2)=  3826  IOR(3)=  7651 

DATAO  After  syntliesls 


DATA: 

.7560E+03 

.0600E+02 

.1200E+02 


2 . 1800E+02 
2.7560E+03 


7.3445E+02 

6.7748E+01 

-2.6972E+01 

-1.5284E+02 


1.8961E+02 

2.0854E+02 

2.1361E+02 

2.2239E+02 


ICOMP=  1 
IZ=  0 
IY=  0 

A(I)  after  anal+sync  (first  six) : 
2.7560E+03  7.3445E+02  1.8961E+02 

IY=  1 

A(I)  after  anal+sync  (first  six): 
2.7560E+03  7.3445E+02  1.8961E+02 

IZ=  1 
IY=  0 

A(I)  after  anal+sync  (first  six): 
2.7560E+03  7.3445E+02  1.8961E+02 

IY=  1 

A(I)  after  anal+sync  (first  six): 
2.7560E+03  7.3445E+02  1.8961E+02 

ICOMP=  2 
IZ=  0 
IY=  0 

A(I)  after  anal+sync  (first  six): 
2.7560E+03  7.3445E+02  1.8961E+02 

IY=  1 

A(I)  after  anal+sync  (first  six): 
2.7560E+03  7.3445E+02  1.8961E+02 

IZ=  1 
IY=  0 

A(I)  after  anal+sync  (first  six): 
2.7560E+03  7.3445E+02  1.8961E+02 

IY=  1 

A(I)  after  anal+sync  (first  six) : 
2.7560E+03  7.3445E+02  1.8961E+02 

ICOMP=  3 
IZ=  0 
IY=  0 

A(I)  after  anal+sync  (first  six): 
2.7560E+03  7.3445E+02  1.8961E+02 

IY=  1 

A(I)  after  anal+sync  (first  six) : 
2.7560E+03  7.3445E+02  1.8961E+02 

IZ=  1 
IY=  0 


A(I)  after  anal+sync  (first  six): 
2.7560E+03  7.3445E+02  1.8961E+02 

IY=  1 


xtest.log 


1.0200E+02  -5.1000E+01 


1.0200E+02  -5.1000E+01 


2.3708E+02 

3.0735E+01 

-5.7078E+01 

-2.7473E+02 


2.3708E+02 


2.3708E+02 


2.3708E+02 


2.3708E+02 


2.3708E+02 


2.3708E+02 


2.3708E+02 


2.3708E+02 


2.3708E+02 


2.3708E+02 


2.3708E+02 


2.0161E+02 

2.1039E+02 

2.1546E+02 

2.3439E+02 


2.0161E+02 


2.0161E+02 


2.0161E+02 


2.0161E+02 


2.0161E+02 


2.0161E+02 


2.0161E+02 


2.0161E+02 


2.0161E+02 


2.0161E+02 


2.0161E+02 


1.0300E+02 


1.0300E+02 


1.2344E+02 

9.3786E-01 

-9.4925E+01 

-8.7585E+02 


1.2344E+02 


1.2344E+02 


1.2344E+02 


1.2344E+02 


1.2344E+02 


1.2344E+02 


1.2344E+02 


1.2344E+02 


1.2344E+02 


1.2344E+02 


1.2344E+02 


^  94/09/19 
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A(I)  after  anal+sync  (first  six): 
2.7560E+03  7.3'J45E+02  1.8961E  +  02 

Just  after  DO  40  loop 
IOR(l)=  1  IOR(2)=  3826  IOR(3)=  7651 

DATAO  just  before  error  chec)<.  loop 


DATA: 

2.7560E+03 

2.0600E+02 

2.1200E+02 

2.1800E+02 

2.7560E+03 


7.3445E+02 

6.7748E+01 

-2.6972E+01 

-1.5284E+02 


1 .8961E+02 
2.0854E+02 
2.1361E+02 
2.2239E+02 


2.3708E+02 


2.3708E+02 

3.0735E+01 

-5.7078E+01 

-2.7473E+02 


2.0161E+02 


2.01G1E+02 

2.1039E+02 

2.1546E+02 

2.3439E+02 


INC1=  1  INC2=  25  INC3=  425 

IOR(l)=  1  IOR(2)=  3826  IOR(3)=  7651 

INDEX=  1  A(INDEX)=  2756.00  DATA(IX)=  2756.00 

INDEX=  3826  A(INDEX)=  2756.00  DATA{1X)=  2756.00 

INDEX=  7651  A(INDEX)=  2756.00  DATA(IX)=  2756.00 


DATA: 

2.7560E+03 

2.0600E+02 

2.1200E+02 

2.1800E+02 

2.7560E+03 


7.3445E+02 

6.7748E+01 

-2.6972E+01 

-1.5284E+02 


1 .8961E+02 
2.0854E+02 
2.1361E+02 
2.2239E+02 


2.3708E+02 

3.0735E+01 

-5.7078E+01 

-2.7473E+02 


2.0161E+02 

2.1039E+02 

2.1546E+02 

2.3439E+02 


FINAL  SUMMARY 


Backtransf orm  equals  original  data 

ICOUNT,  ICOMP,  IX,  lY  and  IZ  at  maximum  error: 

0  0  0  0  0 

CPU  time  for  FFT  routine 

Temperton  FFT77  (anal  +  syn) /2  =  1.333333E-02  secs 


rfesunh%  xtest 
??  DIR=1,  2,  or  3? 

2 

IDIR  =  2 

LBLAS (I)  : 

24  16  8  1  25  425  3825  0  30  0 

Nl=  24  N2=  16  N3=  8 

INC1=  1  INC2=  25  INC3=  425 

IOR(l)=  1  IOR(2)=  3826  IOR(3)=  7651 

DATAO  Values  as  set 


DATA: 


ur\i.  . 

1.8720E+03 

4 .8857E+02 

1.9834E+02 

1.4441E+02 

2.0400E+02 

2.0634E+02 

1.1647E+01 

2.0800E+02 

-2.9726E+01 

2.0966E+02 

2.1200E+02 

-1.6689E+02 

2.1766E+02 

-5.5712E+02 

1.8720E+03 

ICOMP=  1 

IZ=  0 

IY=  0 

Initial 

data  (first  six) 

1.8720E+03 

1.8720E+03 

1.8720E+03 

1 .8720E+03 

1.8720E+03 

IY=  1 

Initial 

data  (first  six) 

4.8857E+02 
IZ=  1 

IY=  0 

4 .8857E+02 

4 .8857E+02 

4.8857E+02 

4.8857E+02 

Initial 

data  (first  six) 

1.8720E+03 

1.8720E+03 

1.8720E+03 

1.8720E+03 

1.8720E+03 

IY=  1 


1 .2344E+02 


1.2344E+02 

9.3786E-01 

9.4925E+01 

8.7585E+02 


1.2344E+02 

9.3786E-01 

9.4925E+01 

8.7585E+02 


6.0044E+01 

7.8937E+01 


1.8720E+03 


4 .8857E+02 


1.8720E+03 


94/0W19 
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Initial  data  (first  six)  : 
4.8857E+02  4.8857E+02  4.8857E+02 

ICOMP=  2 
IZ=  0 
IY=  0 

Initial  data  (first  six)  : 
1.8720E+03  1.8720E+03  1.8720E+03 

IY=  1 


Initial  data  (first  six) : 
4.8857E+02  4.8857E+02  4.8857E+02 

IZ=  1 
IY=  0 

Initial  data  (first  six): 
1.8720E+03  1.8720E+03  1.8720E403 

IY=  1 


Initial  data  (first  six): 
4.8857E+02  4.8857E+02  4.8857E+02 

ICOMP=  3 
IZ=  0 
IY=  0 

Initial  data  (first  six) : 
1.8720E+03  1.8720E403  1.8720E403 

IY=  1 

Initial  data  (first  six)  : 
4.8857E402  4.8857E402  4.8857E402 

IZ=  1 
IY=  0 

Initial  data  (first  six) ; 
1.8720E403  1.8720E403  1.8720E403 

IY=  1 


Initial  d. 
4 .8857E402 
DATAO  Before 

DATA: 

1.8720E403 

2.0634E402 

2.1200E402 


ta  (first  si 
4 .8857E402 
analysis 


4 .8857E402 
1.1647E401 
-1.6689E402 


)  : 

4.8857E402 


1.9834E402 

2.0800E402 

2.1766E402 


IDIR=  2  Just  before  first  CALL  VECFFT 
Just  after  analysis 

IOR(l)=  1  IOR(2)=  3826  IOR(3)=  7651 

IDIR=  2  Just  after  first  CALL  VECFFT 

DATAO  After  analysis 


DATA: 

1.8720E403 

2.0634E402 

2.1200E402 

ICOMP=  1 
IZ=  0 
IY=  0 


4.8857E402 

1.1647E401 

-1.6689E402 


1.9834E402 

2.0800E402 

2.1766E402 


six)  : 

2.0000E402 


six)  : 

1.0100E402 


six)  : 

2 .0000E402 


A(I)  after  anal  (first 
2.0000E402  2.0000E402 

IY=  1 

A(I)  after  anal  (first 
1.0100E402  1.0100E402 

IZ=  1 
IY=  0 

A(I)  after  anal  (first 
2.0000E402  2.0000E402 

IY=  1 


xtestlog 


4.8857E402 


1.8720E403 


4.8857E402 


1.8720E403 


4.8857E402 


1.8720E403 


4.8857E402 


1.8720E403 


4 .8857E402 


1.4441E402 

-2.9726E401 

-5.5712E402 


1.4441E402 

-2.9726E401 

-5.5712E402 


2.0000E402 


1.0100E402 


2.0000E402 


4.8857E402 


1.8720E403 


4.8857E402 


1.8720E403 


4.8857E402 


1.8720E403 


4.8857E402 


1.8720E403 


4.8857E402 


2.0400E402 

2.0966E402 

1.8720E403 


2.0400E402 

2.0966E402 

1.8720E403 


2.0000E402 


1.0100E402 


2.0000E402 


4 .8857E402 


1.8720E403 


4 .8857E402 


1.8720E403 


4.8857E402 


1.8720E403 


4.8857E402 


1.8720E403 


4 .8857E402 


6.0044E401 

-7.8937E401 


6.0044E401 

-7.8937E401 


2.0000E402 


1.0100E402 


2.0000E402 
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A(I)  after 
1 .OlOOE+02 

ICOMP=  2 
12=  0 
IY=  0 

A(I)  after 
2.0000E+02 
IY=  1 

A(I)  after 
1 .OlOOE+02 
IZ=  1 
IY=  0 

Ad)  after 
2.0000E+02 
IY=  1 

Ad)  after 
l.OlOOE+02 

ICOMP=  3 
12=  0 
IY=  0 

Ad)  after 
2.0000E402 
IY=  1 

A(I)  after 
l.OlOOE+02 
IZ=  1 
IY=  0 

Ad)  after 
2.0000E+02 
IY=  1 

Ad)  after 
l.OlOOE+02 


anal  (first 
l.OlOOE+02 


anal  (first 
2.0000E+02 


anal  (first 
l.OlOOE+02 


anal  (first 
2.0000E+02 


anal  (first 
l.OlOOE+02 


anal  (first 
2.0000E+02 


anal  (first 
1  .OlOOE  +  02 


anal  (first 
2.0000E+02 


anal  (first 
l.OlOOE+02 


six)  : 

1 .OlOOE+02 


six)  : 

2.0000E+02 


six)  : 

1 .OlOOE+02 


six)  : 

2.0000E+02 


six)  : 

1 .OlOOE+02 


six)  : 

2 .OOOOE+02 


six)  ; 

l.OlOOE+02 


six)  : 

2.0000E+02 


six)  : 

1 .OlOOE+02 


IDIR=  2  Just  before  second  CALL  VECFFT 
Just  after  synthesis 

IOR(l)=  1  IOR(2)=  3826  IOR(3)=  7651 


DATAO  After  synthesis 


DATA: 

1.8720E+03  4.8857E+02  1.9834E+02 

2.0634E  +  02  l.ie'JYE  +  Ol  2.0800E+02 

2.1200E+02  -1.6689E+02  2.1766E+02 

ICOMP=  1 
12=  0 
IY=  0 

A(I)  after  anal  +  sync  (first  six): 
1.8720E+03  1.8720E+03  1.8720E+03 

IY=  1 


Ad)  after  anal  +  sync  (first  six): 
4.8857E+02  'J.8857E+02  4.8857E+02 

12=  1 
IY=  0 


Ad)  after 
1.8720E+03 
IY=  1 


anal  +  sync  (first  six)  : 
1.8720E+03  1.8720E+03 


A(I)  after 
4.8857E+02 


anal+sync  (first  six): 
4.8857E+02  4.8857E+02 


xtest.log 

l.OlOOE+02  l.OlOOE+02 


2.0000E+02  2.0000E+02 


l.OlOOE+02  l.OlOOE+02 


2.0000E+02  2.0000E+02 


l.OlOOE+02  l.OlOOE+02 


2.0000E+02  2.0000E+02 


l.OlOOE+02  l.OlOOE+02 


2.0000E+02  2.0000E+02 


l.OlOOE+02  l.OlOOE+02 


1.4441E+02  2.0400E+02 
-2.9726E+01  2.0966E+02 
-5.5712E+02  1.8720E+03 


1.8720E+03  1.8720E+03 


4.8857E+02  4.8857E+02 


1.8720E+03  1.8720E+03 


4.8857E+02  4.8857E+02 


l.OlOOE+02 


2.0000E+02 


l.OlOOE+02 


2.0000E+02 


1 .OlOOE+02 


2.0000E+02 


1 .OlOOE+02 


2.0000E+02 


1 .OlOOE+02 


6.0044E+01 

7.8937E+01 


1.8720E+03 


4 .8857E+02 


1.8720E+03 


4 .8857E+02 


ICOMP=  2 
12=  0 
IY=  0 


94/09/19 

09:27:55 


xtest.log 


A(I)  after  anal+sync  (first  six): 

1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03 

1Y=  1 

A(I)  after  anal+sync  (first  six): 

4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02 

IZ=  1 
IY=  0 

A(I)  after  anal+sync  (first  six): 

1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03 

IY=  1 

A(I)  after  anal+sync  (first  six): 

4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02 

ICOMP=  3 
IZ=  0 
IY=  0 

A(I)  after  anal+sync  (first  six): 

1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03 

IY=  1 

A(l)  after  anal+sync  (first  six): 

4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02 

IZ=  1 
IY=  0 

A(I)  after  anal+sync  (first  six): 

1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03  1.8720E+03 

IY=  1 

A(I)  after  anal+sync  (first  six): 

4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02  4.8857E+02 

Just  after  DO  40  loop 
IOR(l)=  1  IOR(2)=  3826  IOR(3)=  7651 

DATAO  just  before  error  check  loop 

DATA: 

1.8720E+03  4.8857E+02  1.9834E+02  1.4441E+02  2.0400E+02  6.0044E+01 

2.0634E+02  1.1647E+01  2.0800E+02  -2.9726E+01  2.0966E+02  -7.8937E+01 

2.1200E+02  -1.6689E+02  2.1766E+02  -5.5712E+02  1.8720E+03 

INC1=  1  INC2=  25  INC3=  425 

IOR(l)=  1  IOR(2)=  3826  IOR(3)=  7651 

INDEX=  1  A(INDEX)=  1872.00  DATA(IX)=  1872.00 

INDEX=  3826  A(INDEX)=  1872.00  DATA(IX)=  1872.00 

INDEX=  7651  A(INDEX)=  1872.00  DATA(IX)=  1872.00 

DATA: 

1.8720E+03  4.8857E+02  1.9834E+02  1.4441E+02  2.0400E+02  6.0044E+01 

2.0634E+02  1.1647E+01  2.0800E+02  -2.9726E+01  2.0966E+02  -7.8937E+01 

2.1200E+02  -1.6689E+02  2.1766E+02  -5.5712E+02  1.8720E+03 


FINAL  SUMMARY 


Backtransform  equals  original  data 

ICOUNT,  ICOMP,  IX,  lY  and  IZ  at  maximum  error: 

0  0  0  0  0 

CPU  time  for  FFT  routine 

Temperton  FFT77  (anal  +  syn) /2  =  1.104167E-02  secs 


rfesunh%  xtest 

??  DIR=1,  2,  or  3? 

3 

IDIR  =  3 

LBLAS (I) : 

24  16  8  1  25  425  3825  0  30  0 


94/09/19 

09:27:55 


xtest.log 


Nl=  24  N2=  16  N3=  8 

INC1=  1  INC2=  25  INC3=  425 

IOR(l)=  1  IOR(2)=  3826  IOR(3)=  7651 

DATAO  Values  as  set 

DATA: 

1.0200E+03  2.3542E+02  2.0200E+02  3.7078E+01 

2.0600E+02  -2.5708E+02  1.0200E+03 


ICOMP=  1 
IZ=  0 
IY=  0 

Initial  data  (first  six): 

1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03 

IY=  1 

Initial  data  (first  six): 

1.0200E+03  1.020CE+03  1.0200E+03  1.0200E+03 

IZ=  1 
IY=  0 

Initial  data  (first  six)  : 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

IY=  1 

Initial  data  (first  six): 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

ICOMP=  2 
IZ=  0 
IY=  0 

Initial  data  (first  six)  : 

1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03 

IY=  1 

Initial  data  (first  six)  : 

1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03 

IZ=  1 
IY=  0 

Initial  data  (first  six): 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

IY=  1 

Initial  data  (first  six)  : 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

I COMP =  3 

IZ=  0 
IY=  0 

Initial  data  (first  six): 

1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03 

IY=  1 

Initial  data  (first  six)  : 

1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03 

IZ=  1 
IY=  0 

Initial  data  (first  six): 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

IY=  1 

Initial  data  (first  six)  : 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

DATAO  Before  analysis 

DATA: 

1.0200E+03  2.3542E+02  2.0200E+02  3.7078E+01 

2.0600E+02  -2.5708E+02  1.0200E+03 


2.0400E+02 


1.0200E+03 


1.0200E+03 


2.3542E+02 


2.3542E+02 


1.0200E+03 


1.0200E+03 


2.3542E+02 


2.3542E+02 


1.0200E+03 


1.0200E+03 


2.3542E+02 


2.3542E+02 


2.0400E+D2 


4 .7421E+01 


1.0200E+03 


1.0200E+03 


2.3542E+02 


2.3542E+02 


1 .0200E+03 


1.0200E+03 


2.3542E+02 


2.3542E+02 


1 .0200E+03 


1.0200E+03 


2.3542E+02 


2.3542E+02 


4.7421E+01 


IDIR=  3  Just  before  first  CALL  VECFFT 


xtest.log 


94/09/19 

iilisiiiiiiii 

Just  after  analysis 

IORa)=  1  I0R(2)=  3826  I0R(3)=  7651 

IDIR=  3  Just  after  first  CALL  VECFFT 

DATAO  After  analysis 

DATA: 


1 .0200E+03 
2.0600E+02  - 

2.3542E+02 

-2.5708E+02 

2.0200E+02 

1.0200E403 

3.7078E+01 

2.0400E+02 

-4 .7421E+01 

ICOMP=  1 

IZ=  0 

IY=  0 

A(I)  after 
2.0000E+02 
IY=  1 

anal  (first 
2.0000E+02 

six)  : 

2.0000E+02 

2.0000E+02 

2.0000E+02 

2.0000E+02 

A(I)  after 
2.0000E+02 
IZ=  1 

IY=  0 

anal  (first 
2.0000E+02 

six)  : 

2.0000E+02 

2.0000E+02 

2.0000E+02 

2.0000E+02 

A(I)  after 
l.OlOOE+02 
IY=  1 

anal  (first 
l.OlOOE+02 

six)  : 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

A{I)  after 
l.OlOOE+02 

anal  (first 
l.OlOOE+02 

six)  : 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

ICOMP=  2 

IZ=  0 

IY=  0 

A(I)  after 
2.0000E+02 

IY=  1 

anal  (first 
2.0000E+02 

six)  : 

2.0000E+02 

2.0000E+02 

2.0000E+02 

2.0000E+02 

A(I)  after 
2.0000E+02 
IZ=  1 

IY=  0 

anal  (first 
2.0000E+02 

six)  : 

2.0000E+02 

2.0000E+02 

2.0000E+02 

2.0000E+02 

A(I)  after 
l.OlOOE+02 
IY=  1 

anal  (first 
l.OlOOE+02 

six)  : 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

A(I)  after 
l.OlOOE+02 

anal  (first 
l.OlOOE+02 

six)  : 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

I COMP =  3 

IZ=  0 

IY=  0 

A(I)  after 
2.0000E+02 
IY=  1 

anal  (first 
2.0000E+02 

six)  : 

2.0000E+02 

2.0000E+02 

2.0000E+02 

2.0000E+02 

A(I)  after 
2.0000E+02 
IZ=  1 

IY=  0 

anal  (first 
2  .OOOOE  +  02 

six)  : 

2.0000E+02 

2.0000E+02 

2.0000E+02 

2.0000E+02 

A(I)  after 
l.OlOOE+02 
IY=  1 

anal  (first 
l.OlOOE+02 

six)  : 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

A(I)  after 
l.OlOOE+02 

anal  (first 
l.OlOOE+02 

six)  : 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

l.OlOOE+02 

IDIR=  3  Just  before  second  CALL  VECFFT 
Just  after  synthesis 

I0R(1)=  1  IOR(2)=  3826  IOR(3)=  7651 

DATAO  After  synthesis 
DATA: 

1.0200E+03  2.3542E+02  2.0200E+02  3.7078E+01  2.0400E+02  -4.7421E+01 


xtest.log 

2.0600E+02  -2.5708E+02  1.0200E+03 

ICOMP=  1 
IZ=  0 
IY=  0 

A(I)  after  anal+sync  (first  six): 

1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03 

IY=  1 

A(I)  after  anal+sync  (first  six): 

1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03 

IZ=  1 
IY=  0 

A(I)  after  anal+sync  (first  six): 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

IY=  1 

Ad)  after  anal  +  sync  (first  six): 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

I COMP =  2 

IZ=  0 
IY=  0 

A(I)  after  anal+sync  (first  six): 

1.0200E+03  1.0200E+03  1.02QOE+03  1.0200E+03  1.0200E+03  1.0200E+03 

IY=  1 

A(I)  after  anal+sync  (first  six): 

1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03 

IZ=  1 
IY=  0 

A(I)  after  anal+sync  (first  six): 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

IY=  1 

A(I)  after  anal+sync  (first  six): 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

ICOMP=  3 
IZ=  0 
1Y=  0 

A(I)  after  anal+sync  (first  six): 

1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03  1.0200E+03 

IY=  1 

A(I)  after  anal+sync  (first  six): 

1.0200E+03  1.0200E+03  1.020CE+03  1.0200E+03  1.0200E+03  1.0200E+03 

IZ=  1 


Ad)  after  anal  +  sync  (first  six): 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

IY=  1 

Ad)  after  anal  +  sync  (first  six)  : 

2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02  2.3542E+02 

Just  after  DO  40  loop 
IOR(l)=  1  IOR(2)=  3826  IOR(3)=  7651 

DATA  ( )  just  before  error  chec)<  loop 

DATA: 

1.0200E+03  2.3542E+02  2.0200E+02  3.7078E+01  2.0400E+02  -4.7421E+01 

2.0600E+02  -2.5708E+02  1.0200E+03 

INC1=  1  INC2=  25  INC3=  425 

IOR(l)=  1  IOR(2)=  3826  IOR(3)=  7651 

INDEX=  1  A(INDEX)=  1020.00  DATA(IX)=  1020.00 

INDEX=  3826  A(INDEX)=  1020.00  DATA(IX)=  1020.00 

INDEX=  7651  A(INDEX)=  1020.00  DATA(1X)=  1020.00 

DATA: 

1.0200E+03  2.3542E+02  2.0200E+02  3.7078E+01  2.0400E+02  -4.7421E+01 


94/09/19 

09:27:55 


xtest.log 


94/09/19 
09}2T:S5 

2.0600E+02  -2.5708E+02  1.0200E+03 

-  FINAL  SUMMARY  - 

Backtransf orm  equals  original  data 

ICOUNT,  ICOMP,  IX,  lY  and  IZ  at  maximum  error: 

0  0  0  0  0 

CPU  time  for  FFT  routine 

Temperton  FFT77  (anal  +  syn) /2  =  1.166667E-02  secs 


test 


94/09/39 

rfesunh%  cd  TEST! 
rfesunh%  xtest 

Test  1:  influences  from  padding  layer 
maximum  error:  0. 

Test  2:  curl  grad  x  =  0 
maximum  error:  5.96046E-08 

Test  3:  curl  (w  x  r)  =  2  w 
maximum  error:  0. 

rfesunh%  cd  .  .  /TEST2 
rfesunh%  xtest 

Test  1:  influences  from  padding  layer 
maximum  error:  0. 

Test  2:  curl  grad  x  =  0 
maximum  error:  5.960‘?6E-08 

Test  3:  curl  (w  x  r)  =  2  w 
maximum  error:  0. 

Test  A',  current  preserved 
maximum  error:  0. 

rfesunh%  cd  ../TESTS 
rfesunh%  xtest 

Test  1:  influences  from  padding  layer 
maximum  error:  0. 

Test  2:  curl  grad  x  =  0 
maximum  error:  1.49012E-08 

Test  3:  curl  (w  x  r)  =  2  w 
maximum  error:  0. 

Test  4:  current  preserved 
maximum  error:  0. 


rfesunh% 


86 
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I  PEGGIE  test  data  examples 


94/08/37 

09j41;3$  mug2.uif 

o_mug2 

'NLRES  L  False  for  NEWRUN,True  for  RESET'  F 
mug2 

CHLABl  replace  these  four  lines 

CHLAB2  by  run  labelling  information 
CHLAB3  which  will  appear  at  the  start 
CHLAB4  of  the  NOUT  channel  output 


' CHLAB5 

A 

**Label  available  to  programmer 

1.1 

/ 

/ 

'NDIARX 

I 

**Channel 

for  diary 

1.2 

/ 

/ 

'NIN 

I 

**Current 

input  channel 

1.2 

/ 

/ 

' NLEDGE 

I 

**Channel 

for  restart  records 

1.2 

/ 

/ 

'NONLIN 

I 

**Channel 

for  input-output 

1.2 

f 

/ 

'NOUT 

I 

**Current 

output  channel 

1.2 

t 

/ 

'NPRINT 

I 

**Channel 

for  printed  output 

1.2 

t 

/ 

'NPUNCH 

I 

**Channel 

for  card  output  (or  equivalent) 

1.2 

r 

/ 

'NREC 

I 

**Current 

record  number 

1.2 

/ 

/ 

'NRUN 

I 

*  *Maximum 

number  of  steps 

1.2 

t 

/ 

'NADUMP 

lA 

**Codes  for  array  dumps 

1.9 

f 

/ 

'NPDUMP 

lA 

**Codes  for  dumping  points 

1.9 

f 

/ 

'NVDUMP 

lA 

**Codes  for  dumping  arrays 

1.9 

/ 

/ 

'NLCHED 

L 

**.TRUE. 

if  class  0  report -head  required 

1.9 

t 

/ 

' NLHEAD 

LA 

**.TRUE. 

if  class  1-9  report-heads  required 

1.9 

t 

/ 

'NLOMTl 

LA 

**Class  1 

subprogram  selector 

1.9 

t 

/ 

'NLOMT2 

LA 

**Class  2 

subprogram  selector 

1.9 

t 

/ 

' NLOMT3 

LA 

**Class  3 

subprogram  selector 

1.9 

t 

/ 

'NLREPT 

L 

**.TRUE. 

if  report  required 

1.9 

t 

/ 

BL  inset  part2  uniforral 
BL  sides  annulu  uniforml 
BL  bottom  disc  uniforml 
BL  handle  ubend  uniforml 

BG  part2  polar_to_rectangular_transition  -6  -10  -16 

SPOLRCT  RADCUR=4 . , RADSTR=5 . , THEMIN=30 . , THEMAX=90 . ,  AXMIN=0 . ,  AXMAX=8 . / 
BG  annulu  polar_with_regular_meshing  -6  -50  -16 

SPOLREG  RADINR=4 . , RADOUT=5 . , THEMAX=300 . , AXMAX=8 . / 

BG  disc  polar_with_regular_meshing  -6  -40  -16 

SPOLREG  RADINR=0 . , RADOUT=4 . , THEMAX=3  60 . ,  AXMAX=1 . / 

BG  ubend  polar_with_regular_meshing  -1  -10  -1 
SPOLREG  RADINR=1 . , RADOUT=2 . , THEMAX=180 . , AXMAX=1 . / 

BP  uniforml  uniform 
SUNIFRM  EPSR=1.,RMUR=1./ 

PP  pcond  pcond 
PA  pcond  sameas  pcond 
BC  inset 

W(0,0)  (l,0.125)=bottom:E(0  .,  0  .)  (0.16666667,  1.0) 

W(0, 0.125)  (l,l)=pcond 
S=sldes  :N 
N=sides :  S 

E(0.4,  0.75)  (0.6,  0.875)=handle:N(0,l)  (1,0) 

E(0.4,0.375)  (0.6,  0.5)=handle:S(l,l)  (0,0) 

E(0.,0.)  (0 . 4 ,  1 . 0)  =pcond 
E(0.6,0.)  (1.0,  1.0)=pcond 
E(0.4,0.)  (0.6,  0.375)=pcond 
E(0.4,0.5)  (0 . 6,  0 . 75)  =pcond 
E(0.4,0.875) (0.6, 1.0)=pcond 
U=pcond 
Down  =pcond 
EN 

BC  handle 
W=pcond 

S(l,l)  (0, 0)  =inset  :E  (0 .4,  0.375)  (0.6,  0.5) 

N{0,1)  (l,0)=inset:E(0.4,  0.75)  (0.6,  0.875) 

E=pcond 

U=pcond 

D=pcond 

EN 

BC  sides 

W(0,0)  (l,0.125)=bottom:E(0. 16666667,  0)  (1,  1) 

W(0, 0.125)  (l,l)=pcond 

S=inset  :N 

N=inset :  S 

E=pcond 

U=pcond 

Down  =pcond 

EN 

BC  bottom 
W=pcond 
S=pcond 
N=pcond 


mug2.uif 

E(0.,  0.)  (0. 16666667,  1.0)=inset:W(0,0)  (1,0.125) 
£(0.16666667,0) ( 1 , 1 ) =sides : K ( 0 , 0 ) (1,0.125) 

U=pcond 
Down  =pcond 
EN 

SF  origin 

OR  handle  inset  sides  bottom 


94/08/17 

09:42:35 


94/08/lg 

tl30.uif 

o_tl30 

'NLRES  L  False  for  NEWRUN,True  for  RESET'  F 
tl30 

CHLABl  replace  these  four  lines 


CHLAB2 

by 

run  labelling  information 

CHLAB3 

which  will  appear  at  the  start 

CHLAB4 

of 

the  NOUT  channel  output 

' CHLAB5 

A 

**Label  available  to  programmer 

1.1 

/ 

'NDIARY 

I 

**Channel  for  diary 

1.2 

/ 

'NIN 

I 

**Current  input  channel 

1.2 

/ 

' NLEDGE 

I 

**Channel  for  restart  records 

1.2 

/ 

'NONLIN 

I 

**Channel  for  input-output 

1.2 

/ 

'NOUT 

I 

**Current  output  channel 

1.2 

/ 

'NPRINT 

I 

**Channel  for  printed  output 

1.2 

/ 

'NPUNCH 

I 

**Channel  for  card  output  (or  equivalent) 

1.2 

/ 

'NREC 

I 

**Current  record  number 

1.2 

/ 

'NRUN 

I 

**Maximum  number  of  steps 

1.2 

/ 

' NADUMP 

lA 

**Codes  for  array  dumps 

1.9 

/ 

'NPDUMP 

lA 

**Codes  for  dumping  points 

1.9 

/ 

' NVDUMP 

lA 

**Codes  for  dumping  arrays 

1.9 

/ 

'NLCHED 

L 

**.TRUE.  if  class  0  report-head  required 

1.9 

/ 

'NLHEAD 

LA 

**.TRUE.  if  class  1-9  report-heads  required 

1.9 

/ 

'NLOMTl 

LA 

**Class  1  subprogram  selector 

1.9 

/ 

'NLOMT2 

LA 

**Class  2  subprogram  selector 

1.9 

/ 

'NLOMT3 

LA 

**Class  3  subprogram  selector 

1.9 

/ 

'NLREPT 

L 

**.TRUE.  if  report  required 

1.9 

/ 

SF  courant  number 
SSF  DTMUL=0.4,NRUN=125, 

NDSPS=1, 2, 

NFSPS=1,2, 

NPS=2, 

NTSPS=1,2, 

CFMTPS='grl'  ,  'grl'  , 

MFLFS=3, 1, 

NFLFS==1,  2, 3, 0,  0,0,  0,0,  0,0,  4, 

NFS=2, 

NPSFS=1,2, 

MDODS=l, 1, 

ND0DS=1, 0, 0, 0,  0,0,0,0,0,0,2, 

NDS=2, 

NPSDS=1,2, 

TIMETS=0.,0.20000000E-10,0.15000000E-11,0.,0.,0.10000000E+10,0.,0.20000000E-10,0.15000000E-11,0.,0., 

NPSTS=1,2, 

NTS=2, 

CNLTIM='Time' ,  'Time'  , 

CNSTIM=' t' , ' t' , 

CNUFRQ='GHz' , 'GHz' , 

CNUTIM=' ns' , ' ns' , 

CNURAT=' /ns' , ' /ns' , 

SLIMDO=0., 0.20000001E+02, 0. , 0.20000001E+02, 

MSUDO=8, 2, 

NCOLDO=l, 1, 

NDO=2, 

NDSDO=l,2, 

NSUDO=l,  2,  3,  4,  6,  7,  8,  9,  0,0,  5, 10, 

CNLDO=' axial  distance  (m)',' axial  distance  (m) ' , 

CNSDO='z  (m)','z  (m)  '  , 

FL1MFL=-10. 0000000, 10.0000000, -10.0000000, 10.0000000,-10.0000000,  10.0000000,  0  .,  20.0000000, 

SCLFL=0., 0.,0.,0.,0., -1.0000000, 

NCOLFL=2, 3, 4,2, 

NFL=4, 

NFSFL=1, 1, 1, 2, 

CFLSPC='int  D.dS','int  B.dS','int  j.dS','int  E.dl', 

CNLFL='total  displacement  current' ,' total  magnetic  flux', 'total  electron  current' ,' radial  voltage', 
CNSFL=' Id'  ,  'phi' , ' le' , 'Vr' , 

CNUFL='A'  ,  'Tm2'  ,  'A'  ,  'V'  , 

NPSU=1, 1, 1,  1,5, 6, 6, 6, 6, 10, 

NSU=10, 

NSUBLK=1, 2, 3, 4, 2, 5,  6,  7, 8,  6, 

NSUD=1, 2, 3,  1,2, 3, 1,2, 3, 1,2, 3, 1,2, 3, 1,2, 3, 1,2, 3, 1,2, 3, 1,2,  3, 1,2, 3, 

NSUDOM=l,  1,1, 1,2, 1,1,1, 1,2, 

NSUDS=1,1,1,1,2,1,1,1,1,2, 

NSUFS=1,  1,1, 1,2, 1,1, 1,1, 2, 

NSUO=0, 0,  0,0, 0,0, 0,0,  0,0,  0,0,  0,3, 0,0, 0,0, 0,0,  0,0,  0,0, 0,0,  0,0,  3, 

NSUPS=1,  1,1,1,2,1,1,1,1,2, 

NSUTS=1, 1, 1,1, 2, 1,1, 1,1, 2, 

NSUX=4, 23,  9,  4, 7,  9,  4,  7,  9, 4, 7,  9,  4, 3,  9,  4, 23,  9,  4, 7,  9,  4, 7,  9,  4,  7,  9,  4, 3,  9, 

NXTPSU=6,  0,0,0,10/ 


I 


O.IOOOOOOOE 


tl30.uif 


94/08/lg 
OgjISjlS 

BG  centre  polar_with_regular_meshing  -5  -24  -10 
SPOLREG  RADINR=0 . , RADOUT=0 . 0005 , THEKAX=3 60 . , AXKAX=0 . 001/ 

BG  outer  polar_with_regular_meshing  -5  -8  -10 
SPOLREG  RADINR=0 . 0005, RADOUT=0 . 001 , THEMAX=120 . , AXMAX=0 . 001/ 
BP  uniforml  uniform 
SUNIFRM  EPSR=1 . , RMUR=1 . / 

BP  uniform2  uniform 
SUNIFRM  / 

BL  bl  centre  uniforml 

BL  bla  sameas  bl 

BL  b2  outer  uniforml 

BL  b3  sameas  b2 

BL  b4  sameas  b2 

BL  b2a  sameas  b3 

BL  b3a  sameas  b2 

BL  b4a  sameas  b3 

PP  pcond  perfect_conductor 

SBCS  / 

PP  reswal  resistive_wall 
SBCS  SURFZ=377 . , STHETA=0 .5  / 

PP  aplfld  applied_f ield 

SBCS  DAPLYA(1)=115.5712,DAPLYA(2)=0.,DAPLYA(3)=0.  / 

PP  axis  polar_axis 
SBCS  / 

PA  pcond  sameas  pcond 
PA  reswal  sameas  reswal 
PA  aplfld  sameas  aplfld 
PA  axis  sameas  axis 
BC  b2 

W=  bl:E  (0,  0) (0.333333,1) 

N=b3:S 
S=b4  :N 
E=pcond 
D=aplf Id 
Up  =  b2a:D 
EN 

BC  b3 

W=bl:E(0. 333333,  0)  (0.666667,  1) 

N=b4 :S 

S=b2:N 

E=pcond 

D=aplf Id 

U=b3a:D 

EN 

BC  b4 

W=bl:E(0. 666667,0)  (l.,l) 

N=b2  :S 
S=b3  :N 
E=pcond 
D=aplf Id 
U=b4a;D 
EN 

BC  bl 

E(0,0)  (0. 333333, l)=b2:K 
E(0. 333333,0)  (0.666667,l)=b3:K 
E(0. 666667,  0)  (l.,l)=b4:W 
D=apl f Id 

U(1,0)  (0,l)=bla:D(l,0)  (0,1) 

S=bl  :N 
N=bl :S 
W=axis 
EN 

BC  b2a 

W=  bla:E  (0,  0) (0.333333, 1) 

N=b3a:S 

S=b4a:N 

E=pcond 

D=b2 :U 

Up  =  reswal 

EN 

BC  b3a 

W=bla:E(0. 333333,0)  (0.666667,  1) 

N=b4a  :S 
S=b2a  :N 
E=pcond 
D=b3  :U 
U=reswal 
EN 


tl30.uif 


94/«8/lg 

0S;J5:I5 

BC  b4a 

W=bla:E(0. 666667,0) 

N=b2a:S 

S=b3a:N 

E=pcond 

D=b4 :U 

U=reswal 

EN 

BC  bla 

E(0,  0)  (0.333333, l)=b2a:W 

E  (0.333333,0)  (0 . 666667, 1) =b3a :W 

E(0. 666667,0)  (l.,l)=b4a:W 

D  (1, 0)  (0, 1)  =bl:U  (1, 0)  (0,  1) 

U=reswal 

S=bla:N 

N=bla:S 

W=axis 

EN 

OR  bl  b2  b3  b4  bla  b2a  b3a  b4a 


uni20,uif 


94/09/26 
13:31  :S8 

o_uni20 

'NLRES  L  False  for  NEWRUN,True  for  RESET'  F 
uni20 

Uniform  B  test 

CHLAB2  by  run  labelling  information 

CHLAB3  which  will  appear  at  the  start 
CHLAB4  of  the  NOUT  channel  output 


'CHLAB5  A  **Label  available  to  programmer  1.1  '  / 
'NDIARY  I  **Channel  for  diary  1.2  '  / 
'NIN  I  **Current  input  channel  1.2  '  / 
' NLEDGE  I  **Channel  for  restart  records  1.2  '  / 
'NONLIN  I  **Channel  for  input-output  1.2  '  / 
'NOUT  I  **Current  output  channel  1.2  '  / 
'NPRINT  I  **Channel  for  printed  output  1.2  '  / 
'NPUNCH  I  **Channel  for  card  output  (or  equivalent)  1.2  '  / 
' NREC  I  **Current  record  number  1.2  '  / 
' NRUN  I  **Maximum  number  of  steps  1.2  '  / 
'NADUMP  lA  **Codes  for  array  dumps  1.9  '  / 
'NPDUMP  lA  **Codes  for  dumping  points  1.9  '  / 
'NVDUMP  lA  **Codes  for  dumping  arrays  1.9  '  / 
'NLCHED  L  **.TRUE.  if  class  0  report-head  required  1.9  '  / 
'NLHEAD  LA  **.TRUE.  if  class  1-9  report-heads  required  1.9  '  / 
'NLOMTl  LA  **Class  1  subprogram  selector  1.9  '  / 
'NLOMT2  LA  **Class  2  subprogram  selector  1.9  '  / 
'NLOMT3  LA  **Class  3  subprogram,  selector  1.9  '  / 
'NLREPT  L  **.TRUE.  if  report  required  1.9  '  / 


SF  courant_number 
SSF  DTMUL=0.4, 

NINIT=1, 

MDPART=1, 

NRUN=50, 

NXPTDD=7, 

BUNI(3)=1./ 

BG  typel  quadrilateral_cylinder  -3  -4  -5 

SQUADRI  XQUAD1=0 . , yQUADl=0 . 4 , XQUAD2=0 . 5, YQUAD2=0 . 5, RXMAX=0 . 3, RZMAX=1 .0/ 

BG  type2  quadrilateral_cylinder  -3  -5 

SQUADRI  XQUAD1  =  0 . 2 ,  YQUAD1  =  0 . 5 , XQUAD2  =  0 . 7, YQUAD2  =  0 . 6, RXMAX=0 . 7 , RZMAX= 1.0/ 

BG  typelr  quadrilateral_cylinder  -3  -4  -5 

SQUADRI  XQUAD1=0.2, YQUAD1=0.5,XQUAD2=0.5, YQUAD2=0.5,RXKAX=0.5,RYKAX=0.1,RZKAX=1.0/ 
BG  type2r  quadrilateral_cylinder  -3  -4  -5 

SQUADRI  XQUAD1=0 . , YQUAD1=0 . 6, XQUAD2=0 . 7 , YQUAD2=0 . 6, RXMAX=0 . 5, RYKAX=0 . 1 , RZKAX= 1.0/ 
BP  uniforml  uniform 
SUNIFRM  EPSR=1.,RMUR=1./ 

BP  uniform2  uniform, 

SUNIFRM  / 

BL  bl  typel  uniforml 

BL  b2  type2  uniforml 

BL  b3  typelr  uniforml 

BL  b4  type2r  uniforml 

BL  bla  sameas  bl 

BL  b2a  sameas  b2 

BL  b3a  sameas  b3 

BL  b4a  sameas  b4 

PP  pcond  perfect_conductor 

PP  reswal  reslstive_wall 

SBCS  SURFZ=377.,STHETA=0.5  / 

PP  aplfld  applied_f ield 

SBCS  DAPLYA(1)=115.5712,DAPLYA{2)=0.,DAPLYA(3)=0.  / 

PA  pcond  sameas  pcond 

PA  reswal  sameas  reswal 

PA  aplfld  sameas  aplfld 

BC  bl 

E=b2:W 

D=bla  :U 

U=bla:D 

S=pcond 

N=b4  :S 

W=pcond 

EN 

BC  b2 
W=bl :E 
N=b3  ;S 
S=pcond 
E=pcond 
D=b2a:U 
U=b2a  :D 
EN 

BC  b3 


uni20.uif 


94/09/26 

W=b4  :E 

S=b2  :N 

N=pcond 

E=pcond 

D=b3a  :U 

U=b3a:D 

EN 

BC  b4 
W=pcond 
S=bl :N 
N=pcond 
E=b3 :W 
D=b4a  :U 
U=b'5a:D 
EN 

BC  bla 

E=b2a:W 

D=bl:U 

U=bl  :D 

S=pcond 

N=b4a:S 

W=pcond 

EN 

BC  b2a 

W=bla;E 

N=b3a:S 

S=pcond 

E=pcond 

D=b2  :U 

U=b2:D 

EN 

BC  b3a 

W=b4a:E 

S=b2a:N 

N=pcond 

E=pcond 

D=b3  :U 

U=b3:D 

EN 

BC  b4a 
W=pcond 
S=bla:N 
N=pcond 
E=b3a:W 
D=b4  :U 
U=b4  :D 
EN 

OR  bl  b2  b3  b4  bla  b2a  b3a  b4a 
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J  MIMDPIC  test  data  examples 


94/09/22 

H;39:0$ 

rfesunh%  ~/3DPIC/MIMDPIC/xmimdpic  test40.uif 
Type  number  of  PROCESSES  required,  NPRES 
must  be  .LE.  number  nodes  =  1 

1 

NPRES  input  as  1 
CNOUT=o_test40pl 
End  of  MODIFY,  MYPRES  =  1 

rfesunh% 


test40.log 


test40.uif 


Mmmi 

14j09;06 

o_test40 

'NLRES  L  False  for  NEWRUN,True  for  RESET'  F 
test'iO 

3-d  particle  lattice  test 
warm  electron  plasma 

2x2x2  blocks  of  4  x  4  x  4  elements 
2x2x2  superparticles  per  element  initially 


' CHLAB5 

A 

**Label  available  to  programmer 

1.1 

/ 

/ 

'NDIARY 

I 

**Channel  for  diary 

1.2 

t 

/ 

'NIN 

I 

♦♦Current  input  channel 

1.2 

t 

/ 

'NLEDGE 

I 

♦♦Channel  for  restart  records 

1.2 

t 

/ 

'NONLIN 

I 

♦♦Channel  for  input-output 

1.2 

/ 

/ 

'NOUT 

I 

♦♦Current  output  channel 

1.2 

t 

/ 

'NPRINT 

I 

♦♦Channel  for  printed  output 

1.2 

t 

/ 

'NPUNCH 

I 

♦♦Channel  for  card  output  (or  equivalent) 

1.2 

! 

/ 

'NREC 

I 

♦♦Current  record  number 

1.2 

/ 

/ 

'NRUN 

I 

♦♦Maximum  number  of  steps 

1.2 

t 

5 

' NADUMP 

lA 

♦♦Codes  for  array  dumps 

1.9 

i 

/ 

'NPDUMP 

lA 

♦♦Codes  for  dumping  points 

1.9 

t 

/ 

' NVDUMP 

lA 

♦♦Codes  for  dumping  arrays 

1.9 

t 

/ 

'NLCHED 

L 

♦♦.TRUE,  if  class  0  report-head  required 

1.9 

t 

/ 

'NLHEAD 

LA 

♦♦.TRUE,  if  class  1-9  report-heads  required 

1.9 

t 

/ 

'NLOMTl 

LA 

♦♦Class  1  subprogram  selector 

1.9 

t 

/ 

'NL0MT2 

LA 

♦♦Class  2  subprogram  selector 

1.9 

t 

12*F, 1*T/ 

'NL0MT3 

LA 

♦♦Class  3  subprogram  selector 

1.9 

t 

/ 

'NLREPT 

L 

♦♦.TRUE,  if  report  required 

1.9 

t 

/ 

'CLIGHT 

R 

♦Speed  of  light  (SI) 

2.1 

0 

/ 

'BAPLYD 

R 

♦applied  B  field  (SI) 

2.2 

f 

0. 

,0 

'BCATR 

RA 

♦♦array  of  boundary  condition  attributes 

2.2 

0 

/ 

'EAPLYD 

R 

♦applied  E  field  (SI) 

2.2 

0 

0, 

.0 

'SPATR 

RA 

♦array  of  particle  species  attributes 

2.2 

0 

2E8/ 

'SURFZ 

R 

♦surface  impedance  in  Z0,s 

2.2 

0 

/ 

'XLENl 

RA 

♦dimension  of  block  type  1 

2.2 

0 

1.2E-4, : 

1.2E 

;-4,  1.2E-4/ 

'XLEN2 

RA 

♦dimension  of  block  type  2 

2.2 

0 

/ 

'NB 

lA 

♦no  of  blocks  in  side  for  NCASE=2 

2.2 

0 

2,2,2/ 

'NCASE 

I 

♦select  device  initialisation  case 

2.2 

0 

2 

'NINIT 

I 

♦select  field  initialisation 

2.2 

0 

4 

'NODIM 

I 

♦dimensionality  of  problem 

2.2 

0 

3 

'NOELl 

lA 

♦elements  in  block  type  1  side 

2.2 

0 

4,4,4  , 

'NOEL2 

lA 

♦elements  in  block  type  2  side 

2.2 

0 

/ 

'NPINIT 

I 

♦select  particle  initialisation 

2.2 

0 

4 

'NSPEC 

1 

♦number  of  particle  species 

2.2 

0 

1 

'COUR 

R 

♦max  courant  number 

3.1 

0 

0. 

,95 

'NBLOCK 

I 

♦number  of  blocks 

4.1 

0 

/ 

'MIMDOP 

I 

♦♦if  =1,  diagnosics  for  MIMD,  =2  for  XPATCH 

5.1 

0 

/ 

'NGMAX 

I 

♦Max  no  of  grid  file  frames 

5.1 

0 

/ 

'NOPSEL 

I 

♦♦Select  Output  Sequence 

5.1 

0 

3 

'NSl 

I 

♦♦output  every  NSl  steps 

5.1 

0 

/ 

'NXPTDD 

I 

♦♦select  0.4  EXPERT  diagnostic  dump 

5.1 

0 

3 

g_t  est  40p 1 

3-d  particle  Lattice  test 
warm  electron  plasma 

2x2x2  blocks  of  4  X  4  X  4  elements 
2x2x2  superparticles  per  element  initially 
Date.  22-09-94 
Times  1 4  s  33  s 1 4 


Refs  test40 
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K  PIC3D  test  data  examples 
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Abstract 

A  benchmark  version  of  the  3D  MIMD  PIC  code  is  described,  and 
measurements  are  reported  for  the  Intel  iPSC/860  and  Paragon  XP/S  140. 
These  show  an  almost  linear  speedup  and  sizeup  from  the  one  processor 
performance  on  the  8  and  64  block  problems.  The  512  block  problem 
shows  a  modest  fall  off  from  linear  performance  scaling  for  more  than  256 
processors,  and  the  4096  block  problem  for  more  than  1366  processors. 
We  conclude  that  the  code  can  make  efficient  use  of  many  hundreds  of 
processors  on  a  large  massively  parallel  computer. 


1  Background 

In  order  to  demonstrate  the  parallel  performance  of  the  new  MIMD  PIC  code, 
a  version  has  been  prepared  as  a  benchmark  program  and  run  on  the  Intel 
Paragon  at  the  Sandia  National  Laboratory,  Albuquerque,  and  on  an  Intel 
iPSC/860.  This  benchmark  version  is  known  as  LPM3  (Local  Particle  Mesh 
benchmark  ^  3)  can  be  used  in  a  standard  way  to  test  the  capabilities  of  the 
code  on  other  computers  available  to  the  USAF  (e.g.  the  IBM  SPl/2  facility  at 
Maui).  It  therefore  also  provides  a  standard  performance  comparison  between 
the  different  computers. 


2  LPM3  Benchmark 

The  earlier  LPM2  benchmark  which  was  provided  to  the  USAF  in  1993  took  a 
typical  two-dimensional  MILO  geometry  to  test  the  parallel  capabilities  of  the 
earlier  code  on  modestly  parallel  computers,  say  up  to  32  or  64  processors.  In 
contrast,  LPM3  is  designed  to  exploit  to  the  maximum  the  capabilities  of  the 
much  larger  Massively  Parallel  Systems  (MPPs)  with  more  than  a  thousand 
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processors  that  are  now  becoming  available  to  the  USAF,  such  as  the  Intel 
Paragon  and  IBM  SP/2.  To  do  this,  the  benchmark  problem  chosen  has  to  be 
three  dimensional  in  order  to  contain  enough  parallelism  to  use  a  large  MPP, 
and  must  also  be  freely  and  easily  scalable  in  size  to  expand  to  the  size  of  any 
MPP.  The  problem  must  also  be  physically  realistic  and  representative  of  the 
problems  for  which  the  code  is  designed.  Our  choice  for  LPM3  was  a  triply 
periodic  three-dimensional  electron  plasma. 

The  plasma  space  is  divided  into  blocks,  each  of  which  contains  512  particles 
representing  the  electrons,  and  64  elements  on  which  the  fields  are  calculated. 
From  the  point  of  view  of  load  balancing  on  the  parallel  computer,  the  block  is 
the  smallest  unit  that  can  be  allocated  amongst  the  processors.  The  problem 
size  is  measured  by  number  of  blocks  Nb,  and  four  problem  sizes  have  been 
used  with  respectively  Nb  =  8,64,512,4096.  These  correspond  respectively  to 
numbers  of  particles  Np  =  4A',327v',2.56A',  2M2  where  A'  =  1024  and  =  A'^. 
The  timestep  is  such  that  about  ten  percent  of  the  particles  leave  each  block 
and  enter  neighbouring  blocks  during  a  timestep.  A  run  of  100  timesteps  is 
chosen  as  the  benchmark  test  because  this  can  be  done  in  a  few  minutes  for 
problem  sizes  and  number  of  processors  of  interest  (timestep  per  second  in  the 
range  0.1  to  10),  and  the  conservation  of  total  number  of  particles  is  used  as  a 
validity  check. 

The  program  also  prints  every  ten  steps  the  number  of  particles  in  block 
one  and  the  number  of  particles  leaving  block  one  in  a  timestep.  Because  the 
particles  in  each  block  are  loaded  with  identical  positions  relative  to  the  block 
edges  and  with  identical  velocities,  the  number  of  particles  in  each  block  should 
remain  at  512  (for  every  particle  leaving  there  is  an  image  particle  entering 
the  block).  This  is  generally  true  up  to  90  steps  or  so,  but  small  differences  in 
rounding  errors  in  the  edge  tests  mean  that  the  number  of  particles  in  each  block 
varies  slightly  from  512  at  the  end  of  100  steps.  This  is  not  considered  an  error 
unless  the  total  number  of  particles  in  the  whole  plasma  changes  (i.e.  particles 
are  being  lost  from  the  system  as  a  whole).  This  has  never  been  observed  with 
the  LPM3  code. 

There  are  two  different  versions  of  the  benchmark  code,  which  are  selected 
by  the  value  of  the  last  input  variable  MXPASW.  The  per-patch  version  sends 
a  separate  message  for  every  patch  in  the  system,  and  there  are  9  patches  for 
every  block.  In  the  per-process  version,  on  the  other  hand,  the  patch  messages 
are  assembled  and  sorted  in  a  buffer  so  that  only  one  message  is  sent  to  every 
other  process  to  which  a  given  process  is  attached.  The  per-process  code  may 
send  10  or  100  times  fewer  messages  than  the  per-patch  version.  The  per-process 
version  should  be  significantly  faster  than  the  per-patch  version  on  computers 
with  a  high  message  startup  time  or  latency.  For  computers  with  low  latency 
there  will  be  little  difference  between  the  versions. 


3  LPM3  Results 

All  the  benchmark  measurements  obtained  are  shown  in  Fig.  1.  The  results 
are  expressed  as  Temporal  performance  in  units  of  timestep  per  second.  This 
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is  a  clear  and  unambiguous  metric  which  expresses  exactly  what  a  program 
user  wishes  to  maximise.  We  have  deliberately  not  used  the  sometimes  popu¬ 
lar  metrics  of  Parallel  Speedup  and  Efficiency  because  these  are  not  absolute 
measures  and  cannot  be  used  correctly  to  compare  the  performance  of  differ¬ 
ent  computers.  Neither  have  we  expressed  the  results  as  Megaflop  per  second 
(Mflop/s),  because  this  would  make  the  curves  lie  closely  on  top  of  each  other, 
and  disguise  the  actual  time  of  computation  from  the  reader.  In  our  units  of 
timestep  per  second  (tstep/s),  calculations  which  take  the  same  time  lie  at  the 
same  height  in  the  graph,  and  faster  calculations  which  take  less  time  lie  higher 
in  the  graph.  These  statements  would  not  necessarily  be  true  if  the  results  were 
expressed  as  speedup  or  in  megaflop  per  second. 

For  each  block  size  there  are  two  pairs  of  curves.  The  pair  closest  to  the 
theoretical  dotted  line  is  that  for  the  Intel  Paragon  running  under  the  Sandia 
SUNMOS  1.4.8  operating  system,  and  the  pair  about  60  percent  lower  and 
slower  are  for  the  iPSC/860  running  under  the  NX  release  3.3.2  operating  sys¬ 
tem.  For  each  computer  the  open  symbols  are  the  measured  values  using  the 
per-process  code  and  the  corresponding  filled  symbols  for  the  per-patch  code. 
The  fact  that  the  curves  for  both  versions  are  almost  the  same  on  both  the  com¬ 
puters,  means  that  message  latency  is  not  a  problem  with  the  Paragon  under 
SUNMOS  or  with  the  iPSC/860.  The  latency  has  separately  been  measured 
using  the  COMMSl  ‘pingpong’  benchmark  to  be  about  80/US  for  both  these 
computers. 

The  dotted  lines  are  theoretical  perfect  linear  scaling  predictions  for  the 
Paragon,  calculated  from  the  one-processor  performance  on  the  eight  block 
problem.  This  scaling  assumes  that  performance  is  proportional  to  the  number 
of  processors,  and  inversely  proportional  to  the  problem  size.  Thus  we  are  saying 
that  using  p-times  as  many  processors  should  ideally  allow  one  to  compute  an 
existing  problem  p-times  faster  (the  speedup),  or  solve  a  problem  p-times  bigger 
in  the  same  time  (the  sizeup).  Alternatively,  of  course,  the  increase  in  number 
of  processors  should  be  able  to  be  used  to  obtain  a  combination  of  speedup  and 
sizeup.  The  stepwise  nature  of  the  measured  performance  curves  arises  due  to 
load  imbalance  when  the  number  of  blocks  is  not  exactly  divisible  by  the  number 
of  processors.  If  we  concentrate  attention  on  the  best  performance  figures  which 
correspond  to  perfect  load  balance,  when  every  processor  is  computing  the  same 
number  of  blocks,  we  find  that  the  measured  performance  is  within  80  percent 
of  the  ideal  linear  scaling  except  for  the  512  block  case  for  greater  than  256 
processors,  and  the  4096  block  case  for  more  than  1366  processors.  For  these 
cases  with  larger  numbers  of  processors,  performance  saturation  is  beginning 
to  be  seen. 

For  the  set  of  cases  shown  in  fig.  1  it  is  evident  that  there  is  not  enough 
parallelism  in  the  problem  to  obtain  more  than  about  10  timestep  per  second. 
The  8  block  problem  can  at  most  use  8  processors,  and  there  is  no  further 
possibility  of  using  more  processors  to  speedup  that  problem.  The  best  we  can 
do  with  the  extra  processors  is  to  solve  bigger  problems  in  about  the  same  time, 
and  this  is  seen  for  the  cases  of  A-i,  —  64,512,4096,  none  of  which  exceed  10 
tstep/s  however  many  proces,sors  are  used. 
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4  Final  Remarks  and  Recommendations 

The  measurements  using  the  LPM3  benchmark  version  of  the  new  3D  PIC 
code  show  almost  linear  speedup  and  sizeup  of  performance  with  the  number 
of  processors  up  to  about  256,  after  which  noticeable  performance  degradation 
from  the  ideal  is  observed.  Further  runs  on  the  Intel  Paragon  are  intended 
for  the  same  problem  but  with  an  eight  times  smaller  block  size.  This  should 
allow  more  processors  to  be  used  to  increase  the  performance  of  each  problem 
size.  After  this  it  is  proposed  that  the  same  tests  be  run  on  the  IBM  SPl  and 
SP2  at  the  USAF  Computer  Center  in  Maui.  This  will  give  a  very  interesting 
comparison  between  the  IBM  and  Intel  systems,  as  well  as  exercising  the  new 
code  on  a  different  manufacturer’s  equipment. 

Further  it  is  recommended  that  the  benchmark  code  be  sanitised  of  pro¬ 
prietory  features  so  that  it  can  be  put  in  the  public  domain  and  become  part 
of  an  internationally  recognised  benchmark  set,  such  as  that  recently  setup  by 
the  PARKBENCH  committee.  If  this  were  done,  the  performance  of  the  code 
would  become  known  on  most  new  computers  as  the  manufacturers  and  others 
run  the  standard  benchmarks,  and  enter  the  results  in  the  publically  available 
data-base,  which  will  be  accessible  via  the  Mosaic  interface  to  the  World-Wide- 
Web. 
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Figure  1:  Temporal  Performance  of  the  LPM3  benchmark  measured  in  units  of 
timestep  per  second,  for  four  problem  sizes  and  up  to  ISfO  processors  on  the 
Sandia  Laboratory’s  Intel  Paragon  XP/S  IfO  (upper  pair  of  curves  which  are 
close  to  the  dotted  lines),  and  up  to  64  processors  on  the  Intel  iPSC/860  (lower 
pair  of  curves).  The  open  symbols  are  measured  values  for  the  per-process  code, 
and  the  filled  symbols  are  the  measured  values  for  the  per-patch  code  (see  text). 
The  dotted  lines  are  theoretical  perfect  linear  scaling  predictions  for  the  Paragon, 
calculated  from  the  one-processor  performance  on  the  eight  block  problem.  This 
scaling  assumes  that  performance  is  proportional  to  the  number  of  proces.sors, 
and  inversely  proportional  to  the  problem  size.  The  stepwise  nature  of  the  per¬ 
formance  curves  arises  due  to  load  imbalance  when  the  number  of  blocks  is  not 
exactly  divisible  by  the  number  of  processors. 
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Abstract 

This  note  outlines  the  revised  system  of  internal  units  used  in  the  kernel 
of  the  3-D  general  geometry  PIC  software. 

It  is  sufficient  to  make  Maxwell’s  equations  dimensionless  in  an  orthogonal 
coordinate  system  as  the  resulting  rules  must  apply  to  the  non-orthogonal  case. 
Let  fundamental  units  be  electric  field  Eq,  velocity  c  and  time  T.  It  is  convenient 
to  introduce  the  auxiUiary  quantity 


L  =  cT. 


Note  T  will  normally  be  set  so  that 


(1) 


T  =  At, 

where  At  is  the  timestep  calculated  as 


(2) 


At  = 


2Co 

Cy/ST 


(3) 


where  Co  <  1  is  a  user  specified  constant  and  F  (defined  in  the  Annex  on 
Dispersion  and  Stability)  contains  geometrical  information.  The  basis  vector 


ei  = 


hence  for  the  orthogonal  discrete  system  where 


(4) 


X  =  hix^,  (5) 

so  that  xi  runs  over  0  to  1  along  a  cell  boundary,  it  follows  that 


ei  =  hiXi  (6) 

where  xi  is  the  unit  vector  in  the  1-direction.  From  its  definition  also 

e^=-ki/hi.  (7) 

There  follows  the  results  that 
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D  =  D(j)x,  =  h(,)D^  .  ^  =  D^e,  (8) 

%) 

=  (9) 

^0) 

Suppose  D(j)  to  have  units  Do,  eqs  (8)  and  (9)  imply  that 

[D,]  =  Do/L  (10) 

and 

[D,]  =  DoL,  (11) 

where  the  square  brackets  denote  the  function  that  returns  the  dimensions  of 
the  enclosed  argument. 

The  relations  between  the  dimensions  of  physical  components  follow  from 
the  continuum  Maxwell  equations  and  are 

Do/T  =  Ho/L  =  Jo,  (12) 

DolL  =  po,  (13) 

Bo/T  =  Eo/L,  (14) 

Do  =  €oEo,  Bo  =  poHo-  (15) 

It  also  follows  that 

Jo  =  (oEo/T,po  =  (oEo/L  (16) 

and 

Bo  =  Eo/c,c^  =  \/{po(o)-  (17) 

The  discrete  equations  are  framed  in  terms  of  d\  Hi  etc.  It  suffices  to  treat  a 
single  component  of  each  vector  field.  Recall  that 

d  =  y/^.  (18) 

Now  yfg  has  dimensions  of  volume,  hence  using  (10)  the  units  of  dJ  are  • 
(Do/T),  ie 


[d^]  =  L^oEo  (19) 

It  follows  that 

[Q]  =  in  =  [Hi]  =  [d']  =  L\oEo  (20) 

since  all  quantities  in  the  difference  equations  representing  Gauss’  Law  and 
Ampere’s  equation  must  have  the  same  dimensions.  Since 


W  After 

3 

b  =  y/gB 

(21) 

it  follows  similarly  that 

m  =  L^Eo/c 

(22) 

and  the  discretised  Faraday’s  Law  implies  also 

[Ei]  =  L'^Eolc. 

(23) 

The  dimensions  of  and  foUow  from  the  constitutive  relations 

r^E-t  [-^l]  1 

^  [di]  eoc’ 

(24) 

II 

II 

O 

(25) 

To  obtain  dimensions  for  4>  and  A  use  the  definition 

1 

o 

1 

II 

(26) 

yielding 

[4>]  =  LEo 

(27) 

and 

[Ai]  =  L^Eolc 

(28) 

Lastly  there  are  the  boundary  conditions  to  treat.  Setting  d’  = 

coC”  implies 

[C-]  =  L'^Eo, 

(29) 

As  might  be  expected  on  general  grounds,  since  (foc)  ^  is  the  impedance  of  free 
space,  it  can  be  shown  that  for  the  surface  impedance  of  walls 

[^]  = 

CqC 

(30) 

The  above  set  of  units  is  given  in  terms  of  Eq.  It  is,  however,  convenient  to  have 
an  expression  for  Eq  in  terms  of  more  fundamental  quantities.  Anticipating  PIC 

code  work  we  take 


Eo 


(31) 


where  mg  is  the  mass  of  an  electron  and  |  e  |  is  the  absolute  value  of  the  charge 
on  an  electron.  The  units  employed  are  summarised  in  Table  1. 
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Table  1:  Dimensionless  units  employed  in  the  code.  T  is  a  shorthand  for  At 
obtained  using  eq  (3)  and  is  given  by  eq  (31). 


Quantity 

Units 

E 

Eo 

B 

Eolc 

Length 

cT 

Time 

T 

Velocity 

c 

^oEo{cTf 

b\Ej 

EocT^ 

EocT 

EocT^ 

a 

EicTf 

Z 

l/(eoc) 

cT 

e’ 

l/{cT) 

_ 

l/(foc) 

foc 

Task  1  Final  Report 


177 


C  Dispersion  and  Stability  Analysis  for  Maxwell’s 
Equations 


AEA/TYKB/31876/TN/8 


DISPERSION  AND  STABILITY 
ANALYSIS  FOR  MAXWELL’S 
EQUATION 

W  Arter 


September  1994 


AEA  Technology 
Culham  Laboratory 
Abingdon 
Oxfordshire 
0X14  3DB 


Document  Control  Number:  AEA/TYKB/31876/TN/8 

Date  of  Issue:  September  1994 

Issue  number:  1 

Authorization 

Name 

Signature 

Position 

Prepared  by 

W  Arter 

U3  .  r, 

Project  Scientist 

Checked  by 

J  W  Eastwood 

'  Project  Manager 

Approved  by 

J  W  Eastwood 

BSSili 

© —  United  Kingdom  Atomic  Energy  Authority,  1994  —  © 


DISPERSION  AND  STABILITY  ANALYSIS  FOR 
MAXWELL’S  EQUATION 

W  Arter 
September  1994 


Abstract 

This  note  derives  the  stability  criterion  used  to  specify  the  timestep 
used  for  integration  of  the  field  equations  in  the  3-D  general  geometry 
PIC  code  PIC3D. 

The  discrete  version  of  Maxwell’s  equations  solved  by  P1C3D  may  be  written 


dtV  =  -e^^'^djGf^d^, 

(1) 

dtd}  = 

(2) 

where  the  overbar  on  the  expressions  containing  and  denotes  that  a 
4-point  average  is  used  for  terms  involving  off-diagonal  G^'^ .  Assuming  that 
G^  and  G^  are  uniform  and  constant,  we  seek  a  solution  to  eqs  (1)  and  (2)  of 
form 


(3) 

^in  _  ^io ^ikjx^ 

(4) 

where 

z  = 

(5) 

The  operators  dt  and  dj 

are  such  that  they  can  be  replaced  as  follows 

dt  {z^  —  1), 

(6) 

Oj  ^  ilij^ 

(7) 

where 

t"  •  L 

A  ,•  =  sm  . 

^  2 

(8) 

The  4-point  average  leads  to  the  appearance  of  a  factor 


cos  ^  cos  ^  =  y/{l  -  Af)(l  -  K])  (9) 

wherever  off-diagonal  G  are  present  in  eqs  (1)  and  (2),  and  G  will  denote  G 
with  this  factor  applied. 

Combining  equations  (1)  to  (8)  gives 
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where 


det(/r^p  -  Mp)  =  0, 


(10) 

(11) 

(12) 


and 

A'^’  =  (13) 

Stability  of  the  discrete  system  requires  |  |<  1  which  may  be  shown  to  be 

equivalent  to 


-4</i<0.  (14) 

No  completely  general  result  is  available  for  the  circumstances  under  which 
solutions  of  eq  (10)  satisfy  (14).  We  treat  the  problem  by  considering  a  number 
of  special  cases. 


Case  (i) 

Under  this  assumption 


where 


H 


m;  =  Q\Q 


e 

p 


Po(o9 


(15) 


(16) 


=  (17) 

is  related  to  the  constant  of  proportionality  in  (15)  and  may  be  taken  as 
unity,  and  in  an  isotropic  medium,  pij  is  the  metric  tensor  (so  gij  contains 
factors  of  the  form  (9)). 

If  Q  has  eigenvalues  since  M  oc  by  a  well-known  result  in  linear 
algebra,  the  eigenvalues  of  M  satisfy 

^(P)  =  (9(p))2^.  (18) 

9 

Using  the  Reduce  computer  algebra  package  [1],  it  may  be  shown  that 


=  Kl{h2h:i  -  9l^)  A  Kl{guh^  -  9l^)  A  Kl{h\h2  -  9I2)  (19) 

+  2K\K2{9\292Z  —  512533)  +  2A^iA3(gi2523  -  513522) 

+  2A'2A'3(5i2513  -  5235ii)- 

If  the  twiddles  in  eq  (19)  are  dropped,  the  stability  criterion  of  Lee,  Palandech 
and  Mittra  [2]  is  found.  However,  the  role  of  the  4-point  averages  is  significant. 
We  use  identities  of  the  type 
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and  write 


A'l  =  sin  y  y  =  ^  >  (20) 

=  ^(1  -  cosfci),  (21) 

1  —  K\  “  ^(1  +  cos^i),  (22) 

Ci  —  cos  ki,  Si  =  sin  fcj.  (23) 


In  2  -  D,  assuming  isotropic  permeability  and  permittivity,  it  may  then 
easily  be  shown  that  equivalent  of  eq  (19)  is  extremal  on  the  boundary  of  the 
Ci  domain,  ie  is  least  when  Ci  =  —1  and  Sj  =  0.  Assuming  that  (19)  also  is 
extremal  on  the  boundary,  ie  at 


Cl  =  C2  =  C3  =  -l,Sl  =  S2  =  ss  =  0,  (24) 

the  stability  criterion 

922933  +  gllg33  +  gllg22^  ^  2  (25) 

is  obtained.  Note  that  the  4-point  averages  have  removed  the  explicit  off- 
diagonal  terms  found  in  Lee,  Palandech  and  Mittra;  however  they  stiU  appear, 
implicitly,  in  g. 


Case  (ii)  No  proportionality  between  Gfj  and 

Attention  now  has  to  focus  on  M^,  since  in  general  M  —  QR  does  not  imply 
^(p)  =  q{p)r(p'l  where  are  the  eigenvalues  of  R.  The  determinantal  equation 
is  readily  formed  by  employing  Reduce.  We  use  the  lack  of  off-diagonal  terms 
above  to  support  the  step,  whereby  the  3  —  D  determinantal  equation  is  made 
tractable  by  restricting  attention  to  the  case  where  Gfj  and  Gfj  are  diagonal. 
It  is  convenient  to  denote  the  scaled  diagonal  elements  of  these  tensors  by  Cj 
and  hi  respectively.  Suppose  once  again  that  the  extremal  /i  is  found  on  the 
boundary  of  c,-  space,  then 


(m*)^  T  +  A  —  0  (26) 

where  n  is  scaled  so  that 


and 


Mg 

c'^AA  ’ 


(27) 


G  =  h^e2  +  ^263  4-  /i3ei  4-  hiCs  4-  hie2  4-  ^261 


(28) 
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A  —  (/i2^3  +  hih^  +  /ii/i2)(c2e3  +  6163  +  e\e2).  (29) 

Analysis  of  eq  (26)  shows  that  inequality  (14)  holds  provided  a  technical 
condition  applies  requiring  that  p  be  real,  which  we  neglect,  and  that 

0  <  -  A  <  8.  (30) 

Introducing 

Ui  —  h2ei,a2  —  hie2,  (31) 

ct3  ^-263,04  (32) 


^5  —  —  ^362,  (33) 

eq  (29)  becomes 

A  =  (03  +  fl5)(®4  "h  ^e)  +  ®i(<^3  "h  ^e)  +  02(^4  +  05)  +  01^2  (34) 

The  expression  24  -  A  is  a  quadratic  in  the  a,  and  has  a  local  minimum  when 
all  a,'  are  equal,  say  to  a,  hence 

24-A  =  12Q-9a^  (35) 

The  expression  (35)  is  positive  provided 

«  <  I-  (36) 

Investigation  of  other  extrema  shows  that  they  give  less  restrictive  criteria  than 
(36). 

Introducing  the  in  vacuo  quantity 

^ij  ~  9ijl y/di  (3/ ) 

relation  (36),  remembering  the  scabng  (27),  implies  that 

cAt^/T  <  (38) 

where 

r  =  rnax{Qf^Q^2iSf2^U^Gu^33->Gf3^Ui^22^^3^^33G^2)':  (39) 

where  Q  contains  only  geometrical  and  relative  permittivity  /  permeability  in¬ 
formation  for  a  general  medium. 

Comments  The  criterion  (38)  is  employed  in  the  code.  In  some  special  cases 
it  is  clearly  unduly  restrictive.  It  has  not  been  rigorously  derived.  However 
it  is  plausible  that  it  ensures  stability  in  general  geometry.  Moreover,  since 
general  geometry  and  anisotropy  of  the  medium  both  correspond  to  significant 
off-diagonal  entries  in  Gij,  it  is  plausible  that  (38)  also  ensures  stability  in 
physically  realisable  anisotropic  media. 
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Abstract 

This  Note  contains  a  specification  of  diagnostics  selection,  the  data 
structure  for  the  diagnostics  and  the  algorithms  for  the  diagnostics  collec¬ 
tion  and  output  in  the  3DPIC  codes. 


1  Introduction 

This  note  describes  the  diagnostics  used  by  the  3DPIC  codes.  It  concentrates 
on  timeseries  diagnostics  which  are  compatible  with  the  TSD  file  format  and 

mpictim[l]. 

Each  process  produces  its  own  set  of  diagnostics  files  which  are  then  merged 
by  a  post  processor  to  produce  diagnostics  files  for  the  whole  simulation.  There 
wiU  usually  be  4  TSD  files  for  each  problem  containing  data  for  points,  lines, 
surfaces  and  volumes  which  correspond  to  the  ‘surfaces’  in  the  TSD  file  format. 
There  will  generally  be  no  TSD  style  ‘audit’  data  items. 


2  Specification  of  Diagnostics 

This  section  describes  the  tagged  input  form  for  specifying  diagnostics. 


PS  = 


(file  format, 


FL  FL  FL  DO  DO  DO 

AW 

su  su  su  su  su  su 


Figure  1:  Plot  Set 
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2.1  Plot  Set  (PS) 

A  Plot  Set  is  specified  by  a  PS  tag  followed  by  a  file  format  name  (’tsd’  for  a 
TSD  file),  the  name  of  a  Field  Set,  the  name  of  a  Domain  Set  and  the  name  of 
a  Time  Set. 

The  tagged  format  is: 

PSuf ile-f ormatuField-SetuDomain-SetuTime-Set 

A  TSD  file  is  produced  for  each  Plot  Set  on  each  process  so  that  the  number 
of  TSD  files  is  the  product  of  the  number  of  processes  and  the  number  of  Plot 
Sets.  A  post  processor  will  combine  the  TSD  files  from  different  processes  so 
that  there  is  only  one  merged  file  for  each  Plot  Set. 

If  a  more  flexible  tagged  file  format  is  used  (e.g.  HDF)  there  need  only  be 
one  file  per  process.  If  message  passing  is  added  there  need  only  be  one  file  for 
all  processes. 

2.1.1  Field  Set  (FS) 

The  Field  Set  specifies  a  list  of  Fields  (FL). 

FSyf ieldset-nameyf ield-nameluf ield-name2ufield-name3u. . . 

2.1.2  Field  Specification 
The  FieLd  specification  format  is: 

FLyf ield-nameyspecif ication-string 

The  specification  string  specifies  which  point  quantities,  line  integrals,  sur¬ 
face  integrals  or  volume  integrals  are  required: 

1.  Point  values  of  scalars  may  be  sampled.  The  scalar  quantities  are  El,  E2, 
E3,  HI,  H2,  H3  (components  of  the  electromagnetic  fields)  nl,  n2,  ... 
(the  numbers  of  particles  of  each  species),  B-H/2,  E-D/2  or  j-E.  These 
quantities  are  represented  by  strings  such  as  ’E.D/2’,  ’j.E/2’,  ’B.H/2’. 

2.  Line  integrals  are  J  E  •  dl  or  J  H  •  dl  and  are  represented  by  the  strings 
’int  E.dl’  or  ’int  H.dl’. 

3.  Surface  integrals  are  J  D  •  dS,  J  B  ■  dS  or  f  j  ■  dS  and  are  represented  by 
the  strings  ’  int  D .  dS  ’ ,  ’  int  B .  dS  ’  or  ’  int  j  .  dS  ’ . 

4.  Volume  integrals  are  f  tl^dV  where  tj}  may  be  one  of  nl,  n2,  . . .,  B  •  H/2, 
E-D/2  or  j  •  E.  The  integrals  are  represented  by  expressions  such  as  ’  int 
nl  dV’,  >int  E.D/2  dV’,  ’int  j.E  dV’  etc. 

Vector  point  quantities,  E,  B,  D,  H  and  j,  are  also  allowed  in  general  but 
are  not  supported  for  TSD  files. 
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2.1.3  Domain  Set  (DS) 

The  Domain  Set  specifies  a  list  of  Domains  (DO).  Each  Domain  Set  belongs  to 
exactly  one  Plot  Set. 


DSudomainset-nameudomain-namelLidomain-name2udomain-name3u . • . 


2.1.4  Domain  (DO) 

The  DOmain  specifies  a  hst  of  SUbdomains  (SU).  Each  Domain  belongs  to 
exactly  one  Domain  Set. 


D0udomain-nameusubdomain-naiaelusubdomain-name2usubdomain-name3u . . . 


2.1.5  Subdomain  (SU) 

The  subdomain  specifies  a  cuboidal  subset  of  a  block.  Each  Subdomain  belongs 
to  exactly  one  Domain. 


SUuSubdomain-nameublock-nameuduOuX 

The  block-name  is  the  name  of  a  block,  o  and  x  are  the  coordinates  of  the 
corners  of  the  subdomain  in  local  normalised  curvilinear  coordinates  and  d  is  a 
triple  of  the  integers  (1,  2,  3,  -1,  -2,  -3)  which  gives  the  directions  for  evaluating 
components  and  integrals.  d(l)  gives  the  direction  for  line  integrals,  d(l)  and 
d(2)  give  the  first  and  second  directions  for  the  vectors  in  a  surface.  d(3)  is  the 
normal  direction  of  a  surface. 

The  o  and  x  points  will  be  converted  to  element  numbers  by  the  prepro¬ 
cessor  and  the  block  name  will  be  converted  to  a  block  number.  All  diagnostic 
quantities  will  be  evaluated  at  element  centres. 


2.1.6  Time  Set  (TS) 

The  Time  Set  specifies  the  times  at  which  output  is  required.  Each  Time  Set 
belongs  to  exactly  one  Plot  Set. 

TSutimeset-nameuStart-timeijend-timeuOUtput-intervalu . . . 

The  start  time,  end  time  and  output  interval  are  given  as  real  numbers. 


3  Data  Structure 

The  following  COMMON  blocks  are  used  to  specify  the  diagnostics: 
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3.1  Plot  Set  Data 
/COMPS/ 

I  MAXPS  maximum  number  of  plot  sets 

I  NPS  **number  of  plot  sets 

lA  NFSPS (MAXPS)  **numbers  of  field  sets  in  each  plot  set 

lA  NDSPS (MAXPS)  **numbers  of  domain  sets  in  each  plot  set 

lA  NTSPS(MAXPS)  **numbers  of  time  sets  in  each  plot  set 

RA  WSPSCL(MAXDOS,MAXFLS)  scalar  workspace 
/CHAPS/ 

AA  CFMTPS*32 (MAXPS)  **file  format  for  each  plot  set 

3.2  Field  Set  Data 
/COMFS/ 

I  MAXFS  maximum  number  of  field  sets 

I  MAXFLS  maximum  number  of  fields  per  field  set 

I  NFS  **number  of  field  sets 

lA  MFLFS (MAXFS)  **number  of  fields  in  each  field  set 

lA  NPSFS(MAXFS)  **plot  sets  of  each  field  set 

3.3  Domain  Set  Data 

/COMDS/ 

I  MAXDS 
I  MAXDOS 
I  NDS 

lA  MOODS (MAXDS) 
lA  NPSDS (MAXDS) 

3.4  Time  Set  Data 
/COMTS/ 

I  MAXIS  maximum  number  of  time  sets 

I  NTS  **number  of  time  sets 

RA  TIMETS(6,MAXTS)  **start,  stop,  output,  average  times,  time  zero,  time  unit 

lA  NSTPTS(6,MAXTS)  **start,  stop,  output,  average  steps,  step  zero 

LA  LCOLTS (MAXIS)  collect  flags 

LA  LOUTTS (MAXIS)  output  flags 

lA  NPSTS (MAXIS)  **plot  sets  of  each  time  set 

RA  TOTLST(MAXTS)  last  output  times 

/CHATS/ 

AA  CNLTIM*64(MAXTS)  **long  names  of  time 

AA  CNSTIM+8 (MAXIS)  **short  names  of  time 
AA  CNUTIM*8 (MAXIS)  **units  of  time 

AA  CNURAT*8 (MAXIS)  **units  of  rate  of  change 
AA  CNUFRQ+8 (MAXIS)  **units  of  frequency 


maximum  number  of  domain  sets 
maximum  number  of  domains  per  domain  set 
♦♦number  of  domain  sets 
♦♦number  of  domains  in  each  domain  set 
♦♦plot  sets  of  each  domain  set 
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3.5  Domain  Data 

/COMDO/ 

I  MAXDD 
I  MAXSUD 
I  NDO 

lA  NDSDO(MAXDO) 
lA  NCOLDO(MAXDO) 

RA  SLIMD0(2,MAXD0) 
RA  SCLD0(2,MAXD0) 
/CHADO/ 

AA  CNLD0*64(MAXD0) 
AA  CNSD0*8(MAXD0) 

3.6  Field  Data 

/COMFL/ 

I  MAXFL 
I  NFL 

RA  SCLFL(2,MAXFL) 
lA  NFSFL (MAXFL) 
lA  NCOLFL (MAXFL) 

RA  FLIMFL(2,MAXFL) 
/CHAFL/ 

AA  CFLSPC*64 (MAXFL) 
AA  CNLFL*64 (MAXFL) 
AA  CNSFL*8 (MAXFL) 

AA  CNUFL*8 (MAXFL) 


maximum  number  of  domains 
maximum  number  of  subdomains  per  domain 
**number  of  domains 
♦♦domain  sets  of  each  domain 
♦♦domain  colours 
♦♦domain  arc  length  limits 
♦♦zero  and  unit  of  arc  length 

♦♦long  names  of  domains 
♦♦short  names  of  domains 


maximum  number  of  fields 
♦♦number  of  fields 
♦♦zero  and  unit  of  field 
♦♦field  sets  of  each  field 
♦♦field  colours 
♦♦field  limits 

♦♦field  specification 
♦♦long  names  of  fields 
♦♦short  names  of  fields 
♦♦units  of  fields 


3.7  Subdomain  Data 


maximum  number  of  subdomains 
♦♦number  of  subdomains 
♦♦subdomain  o  element 
♦♦subdomain  x  element 
♦♦subdomain  directions 
♦♦subdomain  block  number 
♦♦subdomain  domain  number 
♦♦subdomain  domain  set  number 
♦♦subdomain  time  set  number 
♦♦subdomain  field  set  nxmber 
♦♦subdomain  plot  set  number 
lA  LORSUW(MAXFLS,MAXSU)  origin  of  subdomain  workspace  storage 
lA  NDMSUW(4,MAXFLS,MAXSU)  dimensions  of  workspace 
lA  LORSUN(MAXSU)  origin  of  subdomain  node  storage 
lA  NDMSUN(4,MAXSU)  dimensions  of  node  workspace 
lA  NXTPSU(MAXSU)  ♦♦next  primary  subdomain 
lA  NPSU(MAXSU)  ♦♦primary  subdomain  of  subdomain 


I  MAXSU 
I  NSU 

I A  NSU0(3,MAXSU) 
lA  NSUX(3,MAXSU) 
lA  NSUD(3,MAXSU) 
lA  NSUBLK (MAXSU) 
I A  NSUDOM (MAXSU) 
lA  NSUDS (MAXSU) 
lA  NSUTS (MAXSU) 
lA  NSUFS (MAXSU) 
lA  NSUPS (MAXSU) 


JV  J  Brealey,  J  W  Eastwood  and  W  Arter 


1 


4  Algorithms 

4.1  Initialisation 

The  initialisation  routine  DIAINT  is  called  by  PIC3D  from  OUTPUT(l),  It  loops 
through  the  Plots  Sets,  opens  the  diagnostic  files  and  writes  the  headers.  It 
also  clears  the  DIAGWS  array.  This  routine  uses  the  data  structure  in  a  top- 
down  manner. 

The  names  of  TSD  files  have  the  format  run-id- PS- process  .tsd,  where  run- 
id  is  the  run  identification  string  CHREFN,  PS  is  the  plot  set  number  and  process 
is  the  process  number. 

4.2  Data  Collection 

Data  is  collected  by  DIACOL  at  the  subdomain  level. 

In  the  first  call  in  each  timestep  it  determines  which  time  sets  should  have 
data  collected  and  sets  the  LCOLTS  flags  appropriately. 

It  then  loops  through  subdomains  and  checks  LCOLTS  for  the  appropriate 
domain  to  see  if  data  should  be  collected  and  it  checks  if  the  block  containing 
the  subdomain  is  on  the  current  processor. 

If  data  should  be  collected  for  the  subdomain  it  loops  through  the  plots 
in  the  plot  set  and  stores  the  appropriate  data  in  the  DIAGWS  array  using 
LORSUW(JFL,  JSU)  as  the  base  address. 

4.2.1  Low  Level  Collection  Routines 

There  is  one  low  level  diagnostic  routine  for  each  field  quantity.  The  low  level 
routines  operate  on  single  subdomains.  The  routines  have  names  such  as  IJEDV 
for  collecting  J  j  •  EidV.  The  routines  may  return  different  amounts  of  data 
depending  on  the  nature  of  the  subdomain. 

4.3  Data  Output 

The  data  output  routine  DIAOUT  is  called  from  0UTPUT(2). 

The  routine  sets  the  LOUTTS  flags  to  indicate  which  time  sets  need  output. 
It  then  loops  through  plot  sets  in  a  top-down  manner  combining  contributions 
from  subdomains,  scaling,  writing  data  and  clearing  the  DIAGWS  array. 

4.4  Finish 

The  routine  DIAFIN  is  called  from  0UTPUT(2)  to  close  any  files  which  may  be 
open. 

5  Extensions  for  GRID  Line  Plots 

The  diagnostics  have  been  extended  to  produce  Ghost  GRID  files  containing 
plots  of  scalars  along  a  mesh  line.  These  plots  are  selected  by  the  ‘grl’  file 
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format.  The  ‘grl’  files  are  named  in  the  same  way  as  ‘tsd’  files  except  that  the 
suffix  is  ‘ .grl’. 

5.1  Connections  Between  Subdomains 


NXTPSU  =  0 


Figure  2:  Connections  Between  Subdomains 

Connections  between  subdomains  must  be  specified  to  produce  line  plots.  There 
are  two  types  of  connections  end-to-end  connections  and  side-to-side  connec¬ 
tions.  The  subdomains  which  are  connected  end-to-end  are  caDed  primary  sub- 
domains  and  the  other  subdomains  are  called  secondary  subdomains.  The  end- 
to-end  connections  specify  the  Line  along  which  data  is  collected.  The  line  goes 
through  the  o-points  of  each  primary  subdomain  and  goes  in  the  3-direction 
within  each  subdomain.  Each  secondary  subdomain  is  connected  to  exactly 
one  primary  domain  (its  ‘parent’)  and  the  elements  in  the  2  and  3-directions 
must  match  the  corresponding  elements  in  the  parent  primary  domain.  The 
secondary  domains  are  required  so  that  the  flux  through  the  surface,  or  the 
integral  along  a  line,  ‘perpendicular’  to  a  path,  can  be  plotted  against  distance 
along  the  path. 

The  first  primary  subdomain  is  given  by  NSUD0(1,  j).  The  variable  NXTPSU(j) 
is  used  to  point  to  the  next  primary  subdomain.  The  the  last  primary  subdo¬ 
main  and  secondary  subdomains  have  this  pointer  set  to  zero.  The  variable 
NPSU(i)  is  used  to  point  from  a  subdomain  to  its  parent  (a  primary  subdomain 
is  its  own  parent). 

5.2  Plotting  parameters 

The  parameters  NCOLDO  and  NCOLFL  are  used  to  specify  colours  for  each  domain 
and  each  field.  The  parameters  SLIHDO  and  FLIMFL  specify  the  limits  for  plot- 
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ting  arc  length  and  field  value,  SOLDO  and  SCLFL  are  used  for  scaling  the  arc 
length  and  field  values. 

5.3  Workspace  Storage 

Workspace  storage  is  contained  in  the  allocated  array  DIAGWS.  Each  primary 
subdomain  has  storage  for  arc  length  for  nodes  at  address  LORSUN(t)  with  upper 
bounds  for  the  first  four  array  dimensions  given  by  NDMSUN(1  :  4,i).  (Only  the 

3.1  component  of  NDMSUN(1  :  4,i)  will  differ  from  1  for  a  single  line  plot). 
Secondary  subdomains  do  not  have  storage  allocated  for  arc  length  to  nodes 
because  the  values  for  their  parent  subdomains  are  used. 

Field  values  are  stored  in  regions  of  workspace  at  address  LORSUW(i)  with 
dimensions  the  first  four  array  dimensions  given  by  NDMSUW(1  :  4,i).  (Only  the 

3.1  component  of  NDMSUW(1  :  4,i)  wiU  differ  from  1  for  a  single  line  plot).  All 
primary  subdomains  have  storage  allocated  to  them  for  field  values  but  only 
secondary  subdomains  on  the  current  process  have  storage  aUocated. 

5.4  Implementation  Notes 

•  The  tsd  graphics  are  a  special  case  of  the  grl  graphics.  Only  the  final 
output  stages  need  differ. 

•  The  grl  graphics  readily  generalise  to  plots  on  surfaces  and  to  plots  in 
volumes. 

•  Node  arrays  are  one  element  bigger  in  each  non-trivial  direction  than  field 
arrays. 
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Abstract 

The  MPICTIM  code  is  designed  to  examine  time  series  data  produced 
bj'  PIC  simulation  codes.  It  uses  an  OSF/Motif  [1]  Graphical  User  In¬ 
terface  (GUI)  instead  of  the  character  based  interface;  this  OSF/Motif 
interface  to  MPICTIM  is  much  easier  than  a  Command  Line  Interface 
(CLI)  for  interactive  use. 


1  Using  MPICTIM 

This  Section  contains  a  task  orientated  description  of  how  to  use  MPICTIM. 
Section  2  describes  the  functions  of  the  different  components  of  MPICTIM.  The 
MPICTIM  application  uses  OSF/Motif  Graphical  User  Interface  components. 
You  should  refer  to  your  system’s  documentation  for  general  information  on 
using  OSF /Motif.  This  document  concentrates  on  aspects  which  are  specific  to 
MPICTIM. 

1.1  Starting  MPICTIM 

To  start  MPICTIM  start  an  xterm  window  and  issue  the  mpictim  command  at 
the  prompt: 

$  mpictim 

Move  and  resize  the  xterm  window  to  the  lower  right  of  the  screen  as  shown 
in  Figure  1.  Your  workstation  screen  should  now  look  like  the  screen  shown  in 
Figure  1.  You  should  consult  your  system  administrator  if  it  looks  substantially 
different. 

1.2  Quitting  MPICTIM 

To  quit  the  application  activate  the  File  pull  down  menu  from  the  menu  bar 
(Figure  2)  and  select  the  Exit  option. 

1.3  Opening  Data  Files 

To  open  a  time  series  data  file  activate  the  File  pull  down  menu  from  the  menu 
bar  (Figure  2)  and  select  the  Open  option  to  pop  up  the  Open  File  dialogue 
window  (Figure  3).  Use  the  Open  File  dialogue  window  load  a  file  by  selecting 
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Figure  1;  Motif  Time  Series  Analysis  Program  MPICTIM 
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a  file  and  activating  the  OK  button.  Opening  a  time  series  data  file  replaces 
all  time  series  data  held  in  the  code.  If  you  have  sequential  data  from  restart 
calculations  which  occupies  several  files  use  should  use  the  Append  option  from 
the  File  menu  to  append  data  without  deleting  data  which  is  already  loaded. 


1.4  Setting  Processing  Options 

The  data  processing  options  are  set  in  the  Processing  Options  area.  The  options 
take  effect  when  a  data  item  is  selected  and  either  of  the  Read  Audit  Data  or 
Read  Surface  Data,  buttons  are  activated. 

The  period  used  by  the  mean,  rms,  variance  and  spectrum  plots  is  set  by 
typing  text  in  the  Period  text  field  and  activating  the  Apply  button  at  the 
bottom  of  the  main  window. 

The  Differentiate  toggle  button  enables  differentiation  of  the  data  before 
other  processing  is  performed. 

The  Processing  option  menu  (Figure  4)  allows  the  user  to  choose  how  the 
data  is  processed. 

The  Integrate  toggle  button  enables  integration  the  data  after  other  pro¬ 
cessing  has  been  performed. 

1.5  Choosing  the  Data  to  be  Processed 

To  choose  the  data  to  be  processed  select  an  audit  variable  from  the  Audit  Data 
scrolling  list  or  select  select  a  surface  and  a  surface  variable  from  the  Surfaces 
and  Surface  Data  scrolling  lists. 

Activate  the  Read  Audit  Data  button  or  the  Read  Surface  Data  button  to 
process  the  data. 
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/export/home/nick/VIPER2/EHV3RZ2/*,tsd 


Directories 


t/export/home/nick/VIPER2/EHV3RZ2/.  I 


/expor t/home/n i ck/V I PER2/EHV3RZ2/ , , 
/expor t/home/n i ck/V I PER2/EMV3RZ2/ . sb 
/expor t/home/n i ck/V I PER2/EHV3RZ2/D AT A 
/expor t/home/n i ck/V I PER2/EHV3RZ2/0LDD AT A 
/expor t/home/n i ck/V I PER2/EMV3RZ2/REF 


Files 


nl081r25,tsd 

nl081r35,tsd 


best0001«tsd 


test0002*tsd 

test0003*tsd 


Figure  3:  File  Selection  Dialogue 
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Figure  4:  Data.  Menu 


1.6  Producing  Plots 

The  user  may  select  the  colour  used  for  plots  using  the  Colour  option  menu 
(Figure  5).  The  plot  label  is  set  in  the  Label  next  to  the  Plot  button.  The  user 
should  set  the  ranges  for  the  plots  in  the  ranges  fields.  The  Apply  button  at  the 
bottom  of  the  main  window  must  be  activated  to  make  the  values  active.  The 
time  range  is  also  used  as  a  window  when  processing  data.  The  Reset,  Suggest 
and  Data  buttons  at  the  bottom  of  the  main  window  may  be  used  to  display 
current  ranges,  suggested  ranges  and  data  ranges  (and  other  parameters). 

Once  the  parameters  are  set  activating  the  Plot  button  displays  the  plot  on 
the  screen. 

1.7  Annotating  Plots 

The  labels  used  for  annotating  the  plot  are  set  in  the  labels  text  fields.  The 
Apply  button  at  the  bottom  of  the  main  window  must  be  activated  to  make 
the  labels  active.  The  Reset,  Suggest  and  Data  buttons  at  the  bottom  of  the 
main  window  may  be  used  to  display  current  labels,  suggested  labels  and  labels 
from  the  data  file  (and  other  parameters). 

Once  the  labels  are  set  activating  the  Annotate  button  annotates  the  plot 
on  the  screen. 

1.8  Saving  and  Erasing  Plots 

To  save  a  plot  activate  the  Save  button.  The  plot  is  saved  to  a  GHOST  GRID 
file  [2]  with  the  name  given  in  the  text  field  immediately  to  the  right  of  the  save 
button. 

To  erase  a  plot  activate  the  Erase  button. 
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Figure  5:  Colour  Menu 
1.9  Viewing  Parameters 

Information  from  the  time  series  data  files  can  be  displayed  in  the  window  from 
which  MPICTIM  was  started  by  activating  one  of  the  options  from  the  View 
pull  down  menu  on  the  menu  bar  (Figure  6). 

2  Components  of  MPICTIM 

This  Section  describes  the  functions  of  the  different  components  of  MPICTIM. 
Section  1  contains  a  task  orientated  description  of  how  to  use  MPICTIM. 

Figure  1  shows  the  two  base  windows  of  the  MPICTIM  application.  The 
main  control  window  is  top  right  and  the  graphics  display  window  is  top  left. 
The  terminal  window  from  which  MPICTIM  was  started  is  bottom  right. 

The  main  control  area  is  shown  in  greater  detail  in  Figure  7.  The  compo¬ 
nents  of  this  window  are  shown  in  greater  detail  in  Figure  7  and  are  described 
in  the  following  subsections.  The  components  are  described  in  order  starting 
from  the  top  left  and  scanning  from  left  to  right  and  from  top  to  bottom. 

2.1  Menu  Bar 

The  menu  bar  contains  two  pull  down  menus:  the  File  menu  and  the  View 
menu. 

2.1.1  File  Menu 

The  File  menu  contains  three  items: 
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Open  pops  up  a  file  selection  dialogue  to  load  data  from  a  time  series  data  file. 
Loading  the  data  file  replaces  any  data  already  loaded. 

Append  pops  up  a  file  selection  dialogue  to  append  data  from  a  time  series  data  file 
to  data  already  loaded.  The  time  series  data  file  must  contain  the  same 
data  items  and  use  the  same  time  sampling  interval  as  the  files  which  are 
already  loaded. 

Exit  quits  the  application. 


2.1.2  View  Menu 


The  View  menu  is  used  for  displaying  information  in  the  window  from  which 
MPICTIM  was  started.  The  View  menu  contains  four  items: 

Audit  Data  displays  a  list  of  audit  variables.  The  list  gives  the  variable  numbers,  the 
labels,  the  units  used  and  the  long  names. 

Surface  Data  displays  a  list  of  surface  variables.  The  list  gives  the  variable  numbers, 
the  labels,  the  units  used  and  the  long  names. 

surfaces  displays  a  list  of  surfaces.  The  list  gives  the  variable  numbers,  the  labels, 
the  default  plot  colours  and  the  long  names. 


Files  displays  a  list  of  the  files  that  have  been  loaded.  The  list  gives  the  file 
name,  the  run  identification  string,  the  sample  range  in  steps,  the  sample 
range  in  units  of  time,  the  hinted  period  in  steps  and  the  sampling  interval. 


2.2  Plot  Controls 

The  plot  controls  region  contains  controls  used  to  plot  processed  data.  Working 
from  left  to  right  the  controls  in  this  area  are: 
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Figure  7:  Main  Window 
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2.2.1  Plot  Button 

Pressing  the  plot  button  plots  the  processed  data  on  the  screen  using  the  applied 
ranges.  The  plot  is  labelled  with  the  text  from  the  Label  field  (immediately  to 
the  left  of  the  plot  button)  and  is  plotted  in  the  colour  set  in  the  Colour  option 
menu. 

2.2.2  Label  Field 

The  label  field  is  immediately  to  the  left  of  the  plot  button.  The  text  in  this 
field  (up  to  8  characters)  is  used  as  the  label  for  the  current  curve  being  plotted. 
This  field  may  can  be  made  blank  if  no  label  is  required. 

2.2.3  Colour  Option  Menu 

This  control  sets  the  colour  used  to  plot  a.  curve.  The  colour  choices  are 
black, red,  green,  blue,  white,  cyan,  magenta,  yellow.  (White  will  not  be  visible 
because  the  background  is  white). 

2.2.4  Annotate  Button 

This  button  annotates  the  plot  with  the  applied  labels  from  the  Labels  Area. 
Two  title  lines  and  labels  for  the  axes  are  added. 

2.2.5  Erase  Button 

This  button  erases  the  plot. 

2.2.6  Save  Button 

This  button  saves  the  plot  to  a  file  with  the  name  given  in  the  file  name  field 
immediately  to  the  left  and  erase  the  plot. 

2.2.7  File  Name  Field 

The  file  name  field  is  immediately  to  the  left  of  the  Save  button.  It  may  contain 
up  to  16  characters  of  text.  It  specifies  the  name  of  the  GHOST  grid  file  in 
which  the  current  plot  will  be  stored.  The  file  is  placed  in  the  directory  in  which 
MPICTIM  was  started. 

2.3  Processing  Controls 

This  area  contains  controls  which  determine  how  a  data,  item  is  processed  when 
it  is  read. 

2.3.1  Period  Field 

This  field  is  used  for  setting  and  displaying  the  period  used  by  the  mean,  rms, 
variance  and  spectrum  processing  options.  The  period  is  set  when  the  Apply 
button  in  the  Buttons  area  at  the  bottom  of  the  application  is  activated.  The 
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hinted  period  from  the  data  file  is  shown  when  the  Data  button  is  activated, 
the  suggested  period  is  shown  when  the  Suggest  button  is  activated  and  the 
current  value  is  when  the  Reset  button  is  activated. 

2.3.2  Differentiate  Toggle  Button 

When  the  Differentiate  Toggle  button  is  active  the  data  will  be  differentiated 
with  respect  to  time  before  other  processing  takes  place. 

2.3.3  Processing  Option  Menu 

This  option  menu  selects  the  type  of  processing  applied  to  the  data  when  it  is 
processed.  There  are  six  options  available: 

data  is  unchanged 

mean  computed  over  current  period 

rms  (root  mean  squared)  values  computed  over  current  period 

variance  computed  of  current  period 

spectrum  computed  for  time  interval  truncated  to  an  integral  multiple  of  the  current 
period 

strobe  plot  computed  for  time  interval  truncated  to  an  integral  multiple  of  the 
hinted  period 

2.3.4  Integrate  Toggle  Button 

When  the  Integrate  Toggle  button  is  active  the  data  will  be  integrated  with 
respect  to  the  independent  variable  after  other  processing  has  taken  place. 

2.4  Audit  Data  Selection  Area 

This  area  is  used  to  select  an  Audit  Data  item  and  process  it  using  the  current 
options.  The  user  selects  the  variable  in  the  scrolling  list,  the  long  label  is  then 
shown  in  the  selection  text  field  and  the  data  is  processed  by  activating  the 
button  labelled  “Read  Audit  Data’’.  The  processed  data  remains  unchanged 
until  another  Audit  Data  item  or  a  Surface  Data  item  is  read. 

2.5  Surface  Data  Selection  Area 

This  area  is  used  to  select  an  Surface  Data  item  and  process  it  using  the  current 
options.  The  user  selects  the  surface  and  the  variable  in  the  scrolling  lists,  the 
long  labels  are  then  shown  in  the  selection  text  fields  and  the  data  is  processed 
by  activating  the  button  labelled  “Read  Surface  Data”.  The  processed  data, 
remains  unchanged  until  another  a  Surface  Data  item  or  Audit  Data  item  is 
read. 
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2.6  Ranges  Area 

This  area  contains  text  fields  for  viewing  and  setting  various  ranges.  The  values 
in  the  text  fields  are  used  to  set  the  ranges  when  the  Apply  button  in  the 
Buttons  area  at  the  bottom  of  the  window  is  activated.  The  text  fields  are 
set  to  the  current  values  when  the  Reset  button  is  activated,  they  are  set  to 
suggested  values  when  the  Suggest  button  is  activated  and  are  set  to  the  data 
ranges  (after  processing)  when  the  Data,  button  is  activated. 

2.6.1  Time  Range 

These  text  fields  may  show  the  time  range  of  the  processed  data,  the  suggested 
time  range  of  the  processed  data,  the  current  time  range  or  a  range  which 
may  be  applied,  depending  on  which  of  the  Apply,  Reset,  Suggest  or  Data 
buttons  have  been  activated.  This  interval  is  used  as  a  data  window  during 
the  processing  of  data  i.e.  data  outside  the  interval  is  discarded.  This  interval 
is  also  used  as  the  time  interval  for  plots  which  have  time  as  the  independent 
variable. 


2.6.2  Frequency  Range 

These  text  fields  may  show  the  frequency  range  of  the  processed  data,  the 
suggested  frequency  range  of  the  processed  data,  the  current  frequency  range 
or  a  range  which  may  be  applied,  depending  on  which  of  the  Apply,  Reset, 
Suggest  or  Data  buttons  have  been  activated.  This  interval  is  used  as  the 
frequency  interval  for  plots  which  have  frequency  as  the  independent  variable. 
It  is  also  used  as  the  interval  for  the  strobe  plots  in  which  data  from  an  interval 
is  processed  to  produce  one  period  of  the  signal  as  a  function  of  time.  (For 
Strobe  plots  the  time  axis  is  in  units  of  the  hinted  period). 

2.6.3  Variable  Range 

These  text  fields  may  show  the  dependent  variable  range  of  the  processed  data, 
the  suggested  dependent  variable  range  of  the  processed  data,  the  current  de¬ 
pendent  variable  range  or  a  range  which  may  be  applied,  depending  on  which  of 
the  Apply,  Reset,  Suggest  or  Data  buttons  have  been  activated.  This  interval 
is  used  as  the  dependent  variable  interval  for  all  plots. 

2.7  Labels  Region 

This  area  contains  text  fields  for  viewing  and  setting  various  labels  used  for 
annotating  the  plots.  The  values  in  the  text  fields  are  used  to  set  the  labels 
when  the  Apply  button  in  the  Buttons  area  at  the  bottom  of  the  window  is 
activated.  The  text  fields  are  set  to  the  current  values  when  the  Reset  button  is 
activated,  they  are  set  to  suggested  values  when  the  Suggest  button  is  activated 
and  values  from  the  data  file  when  the  Data  button  is  activated. 

Two  lines  of  titles  and  labels  for  the  x  and  y  axes  may  be  set  and  viewed. 
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2.8  Buttons  Area 

This  area  contains  four  buttons  which  control  parameters  used  by  the  applica¬ 
tion: 

Apply  This  button  takes  the  parameters  from  the  Ranges  area,  the  Labels  area, 
the  plot  label  field,  the  Colour  option  menu  and  the  Period  field  and 
makes  the  current. 

Reset  This  button  resets  the  Ranges  area  text  fields,  the  Labels  area  text  fields, 
the  plot  label  field,  the  Colour  option  menu  and  the  Period  text  field  to 
the  current  values. 

Suggest  This  button  sets  the  Ranges  area  text  fields,  the  Labels  area  text  fields, 
the  plot  label  field,  the  Colour  option  menu  and  the  Period  text  field  to 
the  suggested  values. 

Data  This  button  sets  the  Ranges  area  text  fields,  the  Labels  area  text  fields, 
the  plot  label  field,  the  Colour  option  menu  and  the  Period  text  field  to 
the  actual  ranges  of  the  data  and  values  from  the  time  series  data  files. 
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Abstract 

This  note  derives  the  conditions  that  need  to  be  applied  to  Maxwell’s 
equations  in  general  geometry  for  certain  types  of  material  surface.  First 
general  analytic  conditions  are  derived,  then  their  implementation  in  the 
PIC3D  code  is  discussed. 


1  INTRODUCTION 


We  are  solving  Maxwell’s  equation  in  3-D  and  considering  the  conditions  that 
need  to  be  applied  at  a  material  surface.  Let  the  surface  be  defined  by  =  con¬ 
stant,  then  x^  and  x^  are  the  co-ordinates  in  the  surface,  sketched  in  Figure  1. 


Figure  1:  Illustrates  the  geometry  considered  when  deriving  boundary  condi¬ 
tions.  The  vectors  e,  and  62  lie  in  the  surface  which  has  normal  e^. 

First  we  compile  some  useful  vector  identities  relevant  to  such  a  surface.  The 
third  section  derives  the  analytic  form  of  the  boundary  conditions,  and  then 
Section  4  discusses  their  application  in  the  code.  Two  main  types  of  boundary 
are  implemented:  the  resistive  wall  (of  which  the  perfect  conductor  is  a  special 
case)  and  one  that  allows  the  electric  field  to  be  specified  on  a  patch. 
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Boundary  Conditions  for  Maxwell’s  Equations 


2  USEFUL  RESULTS 

The  contravariant  basis  e;  is  related  to  the  covariant  basis  e’  as  follows 


ei  X  02  =  e^^, 

(2.1) 

ei  X  03  =  -e^Vs, 

(2.2) 

02  X  03  =  Vs, 

(2.3) 

1  2  ^3 

0  X  0  —  — 

Vs 

(2.4) 

1  3 

0X0=  — - 

Vs 

(2.5) 

VxV  =  -^ 

Vs 

(2.6) 

where  the 

volume  element 

=  [0]  -02  X  03]. 

(2.7) 

Equations 

(2. 1-2.7)  can  be  used  to  derive  expressions  for  the  tangential  part  of 

a  general  vector  f  with  components  /,,  ie  if 

f  =  fie’  (2.8) 

then  the  cross  product  of  f  with  the  vector  normal  to  the  surface  shown  in  Fig  1 
is 

e^  X  f  =  ——{ei(—f2)  +  e2fi}  (2.9) 

Vs 


or  since  also 


f  =  /'e„  (2.10) 

e^xf  =  (2.11) 

where  by  definition 


9'^=e'-e’.  (2.12) 

There  is  an  important  relation  which  explains  how  it  is  possible  to  define  unam¬ 
biguously  ft,  the  components  of  f  tangential  to  a  surface.  In  Cartesian  geometry 
the  expression  is 

ft  =  (n  X  f)  X  n,  (2.13) 

where  n  is  the  unit  vector  normal  to  the  surface,  and  since  eq  (2.13)  is  a  vector 
relation  it  seems  that  in  general  geometry  the  expression  should  be 


^  (e^  X  f)  X 

I 


(2.14) 
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However,  from  eq  (2.9)  is  clear  that 

fi2  =  (e^  X  f)  X  63  =  eVi  +  eV2 


(2.15) 


which  seems  to  provide  a  more  satisfactory  definition  of  f(  since  clearly  £<2  -ei  = 
/i  etc.  However  from  eqs  (2.9)  and  (2.11)  respectively  also 


—  61(522/1  —  P12/2)  +  62(^11/2  —  P12/1)  (2.16) 

=  61(5""/'  -  g^^f)  +  62(5""/"  -  g^^f)-  (2.17) 


The  equivalence  of  fn  and  ft2  may  be  demonstrated  by  dotting  eq  (2.17)  with 
61  and  62  respectively.  First  however,  it  is  worth  relating  explicitly  and 
9ij  — 


Using  eq  (2.17) 


5"  =  (522533  -  gl3)/g, 

(2.18) 

5^^  =  (511533  -  5?3)/5, 

(2.19) 

5^^  =  (511522  -  9x2)1  g-i 

(2.20) 

9^^  =  -(512533  -  5l3523)/5, 

(2.21) 

5^^  =  (512523  -  9X3922)19, 

(2.22) 

5^^  =  -(511523  -  9X29X3)1 9- 

(2.23) 

5 

(2.24) 

Since  gij  is  the  inverse  of  5®-’,  in  particular  the  (1,  3)  entry  in  their  product 
vanishes,  ie 

gng^^  +  gug^^  +  gi3g^^  =  (2.25) 

so  eq  (2.24)  becomes 

fti  •  6i  =  giif^  +  512/^  +  fi'13/^  (2.26) 


=  /i. 


Similarly 


fsi  •  62  —  /2, 


(2.27) 

(2.28) 


establishing  the  equivalence  of  the  in  the  surface.  Note  that  the  expressions 
are  not  equivalent  in  other  directions,  for 


ft2-6"  =  r-5"V3, 


(2.29) 


fti  •  63  =  /s. 


(2.30) 
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3  ANALYTIC  BOUNDARY  CONDITIONS 


We  take  Maxwell’s  equations  to  be 

D  =  VxH-J,  (3.1) 

V-D  =  /),  (3.2) 

B  =  -V  X  E,  (3.3) 

V-B  =  0,  (3.4) 

and  there  are  also  the  volume  fields  defined  as 

d  =  (3.5) 

b  =  V^B,  (3.6) 

Q  =  x/ffP,  (3.7) 

I  =  y/gj.  (3.8) 


The  boundary  conditions  at  a  wall  implied  by  the  source  equations  (3.2)  and 
(3.4)  are  easily  derived  by  the  standard  shrinking  pill  box  argument  where  n  is 
replaced  by  e^/  |  |.  Let  [f]  here  denote  the  jump  in  a  quantity  /  as  a  wall  is 

crossed,  then  (3.4)  implies 

[B-e3]  =  0,  (3.9) 

which  for  reasonable  geometries  implies  that 

[6^]  =  0  (3.10) 


or 

6^  = 

and  if  B®^‘  is  an  externally  applied  field  with  physical  components  B(i) 


texts 


(3.11) 


since 

I  e"  1= 

Similarly,  equation(3.2)  implies 


[D 


I 


—  Psi 


(3.12) 

(3.13) 


(3.14) 


where  ps  is  the  surface  charge  density,  ie 

[d^]  xz  y/g^ps.  (3.15) 


Next,  consider  the  evolutionary  Maxwell’s  equations,  in  particular  (3.1),  since 
eq(3.3)  is  of  the  same  form.  Integrate  (3.1)  over  a  circuit  C  bounding  a  surface 
S  that  points  in  the  ej  direction,  so  that  the  surface  element  is 
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and  the  line  element 


e'^  X  ei  ,  , 

dl  =  ^ - —AL, 

1  X  ei  1 

(3.17) 

where  c  is  the  height  of  the  surface  (its  extent  parallel  to 
elementary  length.  For  this  circuit 

e^)  and  AL  is  an 

J  3  dS 

(3.18) 

implies 

Now  using  2.9 

H  X  •  ei  =  —  9\2Hi), 

(3.19) 

(3.20) 

and  the  orthogonality  of  ei  and  implies 

1  X  ei  1=  y/g^^gn, 

(3.21) 

hence 

[guH2  -  guHi]  =  -\/gg^^Isi- 

(3.22) 

Similarly  by  considering  a  circuit  bounding  a  surface  pointing  in  the  62  direction 

[g22Hi  -  g\2H2]  =  \/gg^^ls2- 

(3.23) 

Equations  (3.22  and  3.23)  simplify  when  Is  =  0  with  the  consequence  that  the 
equivalent  analysis  for  eq  (3.3)  gives  simply 

[E,]  =  [E2]  =  0  (3.24) 


Equations  (3.22  and  3.23)  are  of  little  use  as  they  stand,  since  the  surface 
current  is  normally  specified  via  the  surface  Ohm’s  law 


[Zh  -  E)  =  0 

(3.25) 

Further  manipulations  are 

begun  by  letting 

=  Zls  -  E. 

(3.26) 

the  geometry  considered 

II 

0 

(3.27) 

CO 

1 

II 

(3.28) 

(oy/g 

(3.29) 

thus 
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for  an  isotropic  medium 

with  uniform  permittivity  cq- 

Now  from  eq  (2.11), 

eq  (3.25)  impbes 

^33^2  ^  ^23^3^ 

(3.30) 

g33fl  ^  gl3f3 

(3.31) 

Thus 

ri  1  .u 

(3.32) 

P=  ^  (d^  ^^^^^) 

(3.33) 

Equations  (3.32)  and  (3.33)  simplify  if  either  (a)  the  surface  co-ordinates  are 
orthogonal  implying  512  =  0,  since  from  eqs.  (2.22),  (2.23)  and  (2.2.5) 


ffl3 

9^^ 

923 

9u  ’ 

g33 

922 

(3.34) 


or  (b)  if  the  3-coordinate  is  normal  to  the  surface  implying  pm  =  §23  =  0,  since 
from  eqs  (2.21)  and  (2.23). 


5^3  =  5^3  =  0.  (3.35) 

Equation  (3.24)  also  cannot  be  employed  as  it  stands  since  the  boundary 
electric  field  is  normally  specified  via  a  potential  difference.  In  the  static  case 

E  =  -V<^  (3.36) 

and  if  there  are  no  charges  on  the  boundary,  which  is  assumed  to  have  uniform 
isotropic  permittivity  coi  then 


0 

11 

1> 

(3.37) 

Equation  (3.36)  implies 

E  - 

'  dx}  ’ 

(3.38) 

which  in  eq  (3.37)  yields 

(3.39) 

for  some  constants  C',  or  on 

multiplying  (3.39)  by  g^j 

p  _ 

^3  ~  f— 

y/g 

(3.40) 

Hence  from  the  definition 

P  _  9ijd^ 

(3.41) 

eq  (3.40)  implies 

d'  =  eoC', 

(3.42) 
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where  the  C'  follow  from  the  specified  potential  difference,  as  will  now  be  illus¬ 
trated. 

Suppose,  for  example,  there  is  a  potential  difference  only  in  the  1-direction, 
then  use 

^  (3.43) 

dx^  y/g 

to  obtain 


Ja 


(3.44) 


where  A  and  B  are  points  on  opposite  boundaries  of  a  patch.  Further  assuming 
no  normal  field  =  0)  and  surface  orthogonal  co-ordinates  (pi2  =  0)  gives 
=  0  and 


fB 

JA  y/g 


(3.45) 


4  IMPLEMENTATION  OF  BOUNDARY  CONDI¬ 
TIONS 

The  conditions  of  applied  magnetic  and  electric  field  are  straightforward  given 
that  in  the  latter  case  the  user  is  able  to  calculate  the  C,  from  knowledge  of  the 
potential  difference,  when  the  condition  (3.42)  can  be  directly  enforced.  The 
main  difficulties  arise  in  the  case  of  resistive  walls. 

The  first  point  to  establish  is  the  relation  between  P  employed  in  the  PIC3D 
code  and  /'.  Integrating  over  an  elementary  volume  spanned  by  the  vectors  ei, 
02  and  03  with  01  and  02  taken  to  define  the  surface,  yields 


I  =  I  01  X  02  I 


(4.1) 


=  (4.2) 

(and  similarly  for  the  charge  densities  Q  and  ps).  Thus  from  eqs  (3.32)  and 
(3.33)  the  contribution  to  I  from  the  resistive  wall  is 

Ik  =  Y^(d'=  -  *=1,2,  (4.3) 

assuming  (temporarily)  pointwise  constitutive  relationships.  Equation  (4.3)  can 
be  easily  expressed  in  terms  of  nodal  values  in  either  of  cases  (a)  or  (b):  the 
second  is  trivial,  and  in  case  (a),  g^^d^/g^^  can  be  replaced  by 


^fk)(k)  ’ 


(4.4) 


where  the  overbar  denotes  that  a  four-point  average  is  used  to  evaluate  the 
numerator.  In  eq  (4.4)  the  coefficients  are  needed  at  points  conforming  to  the 
standard  pattern  for  Gfj  values,  so  no  special  Gfj  calculations  are  needed. 
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The  PIC3D  implementation  of  the  resistive  wall  update  involves  the  con¬ 
struction 


=  +  +  1,2  (4.5) 

where  n,  n-f  1  denote  time  levels,  and  d^  denotes  a  contribution  from  V  x  H  to 
the  field  advance.  Using  a  temporally  decentred  representation  for  d^'  in  (4.3), 
viz  introducting  a  parameter  6  so  that 


^ 

(4.6) 

and  recalling  a  factor  of  2At,  it  follows  that 

Z  +  2e-2 

a  —  — ^ - 

Z  +  20 

(4.7) 

Z  +  20 

(4.8) 

and 

25;, 

^  ZE2e 

(4.9) 

where 

'7  _ 

(4.10) 

and  Sk  is  the  discrete  approximation  to  Observe  that  provided 

0  >  1/2  the  above  can  be  used  to  represent  a  perfectly  conducting  wall  {Z  =  0). 
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Abstract 

Reasons  are  put  forward  for  preferring  the  type  of  interface  that  we 
have  selected  for  our  mesh  generation  code.  The  concepts  underlying  our 
approach  to  mesh  generation  are  outlined,  and  then  a  detailed  specification 
is  presented.  This  document  supersedes  an  earlier  note  issued  on  this  topic 

[3] 

1  Introduction 

The  complexity  of  interfaces  to  existing  mesh  generation  software  has  led  us  to 
produce  a  simpler  interface  tailored  to  the  needs  of  the  program  under  develop¬ 
ment.  The  concept  underlying  our  approach  is  that  the  device  to  be  modelled  is 
made  by  joining  together  ‘parts’,  with  appropriate  coding  to  generate  the  mesh 
and  associated  information  for  each  class  of  part.  One  important  advantage 
is  that  ‘appropriate’  coding  could  just  be  that  which  converts  the  output  of 
another  mesh  package  into  a  format  consistent  with  our  new  software. 

For  reasons  of  simplicity  and  portability,  input  to  our  mesh  generator  takes 
the  form  of  a  file  containing  text.  Previous  experience  highlights  the  value  of 
using  tags  within  the  file,  that  take  the  form  of  two  character  words  appearing 
at  the  start  of  an  item.  The  tags  serve  to  teU  the  processor  the  format  of 
subsequent  data.  The  following  discussion  indicates  the  formats  that  may  be 
required. 

First,  let  us  observe  that  each  class  of  part  must  have  a  unique  name  since 
it  may  correspond  to  a  piece  of  code,  which  good  programming  practice  implies 
must  be  a  subroutine.  In  this  case  there  will  be  parameters  to  be  set,  e.g. 
the  dimensions  of  a  block  and  the  number  of  mesh-points.  Several  different 
parts  may  belong  to  the  same  class,  hence  some  sort  of  marking  is  needed  to 
distinguish  one  part  from  another.  Moreover,  there  are  effectively  two  sorts  of 
class,  not  only  classes  which  correspond  to  subprograms  of  the  mesh  software, 
but  also  those  that  are  determined  dynamically  e.g.  by  reading  the  output  of 
another  mesh  generator.  Such  output  may  be  extremely  verbose,  hence  it  is 
desirable  to  be  able  to  teU  the  code  to  look  in  a  disk  file  for  the  data. 

The  characteristics  required  of  the  most  general  tagged  record  are  thus 
that  it  should  contain  a  label,  a  class  name,  some  way  of  specifying  whether 
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a  disk  file  is  to  be  searched,  and  then  some  parameters.  Experience  again 
suggests  that  all  but  very  small  sets  of  parameters  are  best  specified  using  the 
NAMELIST  format  of  input.  NAMELIST  formats  are  defined  as  part  of  the 
FORTRAN  90  standard  [1],  and  most  FORTRAN  77  compilers  of  which  we  are 
aware  accept  NAMELIST,  with  slight  variations  that  can  be  easily  corrected  for. 
NAMELIST  is  not  so  useful  if  aU  parameters  in  a  set  have  to  be  specified,  hence 
both  NAMELIST  and  ordinary  lists  will  be  permitted.  The  resulting  record 
structure  enables  compactness  of  representation  coupled  with  great  flexibility 
and  readabihty. 

The  key  concepts  underlying  our  approach  to  mesh  generation  are  further 
expounded  in  Section  2.  Section  3  describes  the  three  basic  input  formats 
required,  of  which  the  first  or  standard  is  by  far  the  most  important.  Section  4 
contains  a  detailed  specification  for  input  to  the  mesh  generation  software,  and 
the  final  section  5  describes  the  classes  available. 


2  Key  Concepts 

As  already  mentioned,  a  device  is  conceived  of  being  made  up  of  parts.  Let  us 
think  of  a  part  as  relating  to  a  single  multiblock  or  uniblock.  For  book-keeping 
purposes  it  is  necessary  for  every  uniblock  to  have  a  unique  label.  However, 
just  as  say  many  cars  have  four  identical  wheels,  so  several  uniblocks  may 
correspond  to  the  same  part.  Defining  a  uniblock  thus  involves  (i)  giving  the 
uniblock  a  label  and  (ii)  specifying  the  part  to  which  it  relates,  by  giving  the 
part  description. 

The  part  description  is  composed  of  two  labels,  since  it  is  useful  to  be  able 
to  separate  the  part  geometry  from  the  part  physics,  i.e.  to  be  able  to  specify 
permittivities  and  conductivities  etc.  independently  of  the  shape  of  the  part. 
The  part  physics  is  specified  either  by  giving  a  class  and  associated  parameters, 
or  as  an  entity  on  disk.  The  part  geometry  may  be  specified  similarly. 

Having  specified  the  parts,  or  more  precisely  the  multiblocks,  it  is  necessary 
to  describe  how  they  are  put  together.  If  two  blocks  A  and  B  have  a  common 
face,  then  it  is  only  necessary  to  give  the  information  that  a  particular  face  of 
block  A  maps  to  a  face  in  block  B.  The  convention  for  labelling  faces  is  indicated 
in  Fig  1  which  shows  how  the  labels  N,  S,  E,  W,  U  and  D  (respectively  North, 
South,  East,  West,  Up  and  Down)  relate  to  the  coordinates  used  within  a  block 

If  blocks  only  meet  over  part  of  a  face  then  it  is  necessary  to  specify  which 
fraction  of  each  block’s  face  is  involved.  Incorporating  ideas  from  ref.[2],  con¬ 
nectivity  is  specified  via  surface  patches.  Further,  it  is  helpful  to  consider  the 
co-ordinate  system  in  which  the  block  maps  to  the  unit  cube.  (Such  a  system  is 
obtained  by  a  straightforward  contraction  from  the  system  proposed  in  ref. [2]). 
The  unit  cube  system  is  then  used  to  measure  the  extent  of  the  surface  patch. 
The  surface  patch  is  specified  by  giving  the  face  label  and  co-ordinates  of  a 
pair  of  opposite  corners  (the  ‘O’  and  ‘X’  points).  Note  the  implication  that  the 
patch  shape  in  physical  space  is  effectively  determined  by  the  surface  in  which 
it  lies. 

To  specify  connectivity  in  the  input  file,  a  list  rather  than  a  NAMELIST 
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Figure  1;  fllustrates  the  labelling  of  faces  of  a  multiblock.  Position  xq  is  the 
origin  of  the  block  co-ordinates  and  vectors  xj,  X2  and  X3  indicate  the  1,  2  and 
3  co-ordinate  directions  respectively. 


is  appropriate,  since  we  adopt  the  convention  that  for  each  block,  connectivity 
information  for  every  face  must  appear.  Although  this  implies  a  factor  of  two 
redundancy  in  the  input,  since  clearly  the  same  patch  information  appears  in 
the  two  entries  for  the  adjacent  blocks,  validation  of  the  input  is  considerably 
simplified. 

It  is  desirable  also  to  specify  the  order  in  which  the  parts  are  assembled. 
The  main  reason  is  the  need  to  determine  how  many  elements  are  used  in  each 
surface  patch.  It  is  easy  to  see  that  conflicts  can  occur.  For  example  in  Block 
A  there  may  have  to  be  say  12  X  12  element  faces  in  the  face  common  to  Block 
B,  because  Block  A  has  been  gridded  using  a  foreign  mesh  generator,  but  the 
chosen  gridding  for  Block  B  might  be  16  xl6  X  16  elements.  The  order  of 
assembly  can  be  used  to  resolve  conflicts  by  preferring  data  associated  with  the 
block  that  appears  earlier  in  the  ranking.  In  the  example  shown,  if  A  appears 
before  B,  then  the  code  checks  to  see  whether  a  12  X  12  x  12  gridding  is  possible 
for  Block  B,  which  will  usually  be  the  case  if  B  is  specified  as  a  subprogram 
class.  If  B  appears  before  A  then  the  conflict  cannot  be  resolved  and  an  error 
message  should  result. 

The  above  example  serves  to  introduce  the  concept  of  elastic  grid  sizes 
n\  X  712  y.  ^3  where  n,-  is  the  number  of  elements  in  a  coordinate  direction  i. 
For  most  blocks,  the  associated  n,-  are  parameters  that  can  be  changed  to  try 
to  resolve  conflicts.  The  convention  adopted  is  that  instead  of  using  positive 
integers  for  the  n,-  in  these  blocks,  a  negative  integer  is  used  for  each  co-ordinate 
direction  in  which  the  grid  size  is  effectively  arbitrary.  It  is  worth  noting  that 
even  when  all  blocks  have  been  collected  some  n,-  may  remain  elastic,  and 
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means  is  provided  for  defining  these.  The  default  is  to  replace  negative  iii  by 
their  absolute  values,  so  that  in  effect  negative  n,  suggest  grid  sizes  but  do  not 
prescribe  them. 

3  Basic  Input  Formats 

The  Section  begins  by  describing  the  standard  tagged  input  format  then  goes 
on  to  outline  other,  simpler  tagged  formats.  Detailed  specifications  for  each  tag 
are  given  in  Section  4. 

3.1  Standard  Format 

The  standard  format  of  a  tagged  record  as  used  in  Section  4  is  written: 

QQulabeluChowstored=]classuCparlupar2. . .]. 

Labels  are  case  sensitive  and  should  not  contain  spaces,  commas  or  colons. 
If  the  classname  contains  spaces,  these  may  be  replaced  by  underscore  or  the 
entire  name  should  be  enclosed  in  quotes. 

To  explain  the  above,  note  first  the  convention  that  items  appearing  within 
square  brackets  are  optional.  The  three  compulsory  items,  each  separated  by 
at  least  one  space  “u  ”,  are  QQ  which  is  the  generic  form  for  a  tag,  a  label 
which  does  not  occur  in  any  other  record  with  the  same  tag.  and  a  classname. 
An  example  of  such  a  simple  record  is: 

BGupartlguregular_cubic_lattice 

This  creates  a  block  geometry  entry  that  is  referenced  by  the  label  partlg  and 
described  as  belonging  to  the  class  regular_cubic_lattice  so  that  e.g.  there 
exists  a  subroutine  CUBREG  to  generate  the  particular  geometry. 

Parameters  that  determine  parffgr  within  the  class  regular_cubic_lattice 
are  specified  in  subsequent  lines  of  input,  using  NAMELIST  as  defined  in  the 
standard[l].  Thus  the  complete  entry  for  partlg  also  includes,  e.g. 

&CUBREG 

RXMAX  =  0.1,  RYMAX  =  2.0,  RZMAX  =  3.0/ 

that  specifies  the  dimensions  of  the  cuboid  region  (in  metres). 

The  optional  entries  pari,  par2,  ...  are  special  parameters  that  can  be  set 
without  recourse  to  the  verbosity  of  NAMELIST.  Their  meaning  depends  on 
the  tag.  For  BG  they  are  the  n,,  in  order  i  =  1,2,3,  defining  the  grid  size.  Thus 
a  more  complicated  entry  might  be 

BGupartluregular_cubic_latticeu-20|j-60u-40 
The  use  of  minus  denotes  the  size  is  elastic  (see  Section  2),  but  that  a  suggested 
mesh-size  for  the  part  is  20  X  60  X  40  elements. 

As  its  name  suggests  the  optional  howstored  entry  indicates  how  the  infor¬ 
mation  is  stored.  If  it  is  omitted  a  namelist  normally  follows  as  in  the  example. 
However  other  values  are  possible,  depending  on  the  class  and  tag.  AH  tags 
allow  howstored  to  take  the  value  file,  indicating  the  file  with  the  name  class 
is  to  be  searched  for  the  remaining  information  associated  with  the  entry.  The 
file  class  will  itself  contain  an  entry  of  the  above  form. 
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Values  for  howstored  that  may  be  allowed  depending  on  the  tag  are  number 
and  metric.  Number  implies  that  the  block  geometry  is  specified  as  a  set  of 
numbers  representing  e.g.  the  co-ordinates  of  mesh-points.  Metric  indicates 
that  the  metric  information  needed  for  the  timestepping  calculation  is  already 
available.  The  above  serve  directly  or  indirectly  to  define  new  classes  of  geom¬ 
etry. 

3.2  List  of  Labels 

There  are  basically  two  simpler  formats  for  input.  One  consists  of  a  tag  followed 
by  a  list  of  labels,  constituting  a  complete  entry.  This  is  used  e.g.  to  define 
blocks  in  terms  of  predefined  block  geometries  and  physics  thus 
BLybl lupart Igupart Ip 

and 

BLubl2usameasubll . 

Note  that  sameas  has  a  special  meaning  and  should  not  therefore  be  used  as 
a  block  label.  Any  blocks  defined  without  use  of  a  sameas  can  also  be  thought 
of  as  “block  types”.  The  same  format  is  used  to  specify  the  order  of  assembly 
of  blocks. 

3.3  Connectivity  Format 

The  second  simpler  input  format  is  used  to  define  connectivities.  The  tagged 
line  contains  only  a  label  (in  addition  to  the  tag).  Subsequent  lines  contain  the 
block  connectivity  information  as  series  of  entries  (one  to  a  line)  in  the  format: 

facel[{x\^,x^){x\^,x\^)]  =  block  :  face2[{xl'^,xl'^){x^^,xl^)] 

terminated  by  a  line  containing  only  the  tag  EN.  Omitting  the  quantities  in 
square  brackets  is  equivalent  to  taking 

{xf,xf){x\^,x^2^)  =  (0,0)(1,  l),i=  1,2. 

The  and  are  the  co-ordinates  of  the  ‘0’  and  the  ‘X’  points  defining 
a  patch  on  the  face  facel  where  i  =  1  refers  to  the  block  of  which  the  label 
appears  on  the  tagged  line  and  i  =  2  refers  to  other  blocks,  the  labels  of  which 
appear  immediately  to  the  right  of  the  equals  sign. 

4  Detailed  Input  Specification 

The  allowed  tags  are  assigned  one  of  three  input  formats.  The  tags  BG,  BP 
and  PP  are  used  to  define  respectively  the  geometry  of  a  block,  the  physics  of 
a  block,  and  physical  data  associated  with  a  surface  patch  (including  boundary 
conditions).  The  tag  DE  is  used  to  set  default  values  for  parameters,  and  SF 
is  used  to  set  parameters  that  otherwise  have  no  associated  tag.  The  input 
associated  with  all  these  tags  has  the  standard  format. 

The  tags  BL  and  PA  are  used  to  define  respectively  a  block  in  terms  of 
its  geometrical  and  physical  attributes  and  a  patch  in  terms  of  its  physical 
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attributes.  Tag  OR  is  used  to  specify  the  order  in  which  blocks  should  be 
connected.  The  input  associated  with  these  tags  has  the  list  of  labels  format. 

Lastly  the  tag  BC  is  used  to  define  block  connectivity  information.  The 
information  in  Section  3.3  constitutes  a  detailed  input  specification  for  this  tag, 
thus  no  further  description  appears  here. 

The  Table  provides  a  summary  of  formats  for  the  different  tags. 


Table  1 :  Tag  Table 


Tag 

Brief  Description 

Format -f 

Special  Parameter 

Notes* 

BC 

Block  Connectivity 

C 

- 

- 

BG 

Block  Geometry 

S 

iVi,A2,A3 

f  n  m 

BL 

BLock  definition 

L 

- 

3  labels 

BP 

Block  Physics 

S 

NuN2,N3 

f  n  m 

DE 

Defaults  (namelists) 

S 

- 

f 

OR 

ORdering  of  blocks 

L 

- 

any  number  of  labels 

PA 

PAtch  definition 

L 

- 

3  labels 

PP 

Patch  Physics 

S 

NuN2,N3 

f  n 

SF 

Set  Free  parameters 

S 

- 

f 

*f,  n  and  m  are  abbreviations  for  the  allowable  howstored  options  file,  numbers 
and  metric. 

+S,  L  and  C  denote  the  standard,  the  list  of  labels  and  the  connectivity  formats 
respectively 


4.1  Tags  in  Standard  Format 

4.1.1  BG,  BP  and  PP 

Input  associated  with  the  tag  BG  has  been  specified  in  Section  3.1.  Much 
the  same  detailed  specification  applies  to  tag  BP.  The  only  major  differences 
is  that  the  metric  information  now  relates  to  permittivity,  permeability  and 
conductivity,  rather  than  to  position.  As  might  be  expected,  the  associated 
namelist  names  also  differ.  Tag  PP  input  is  very  similar,  with  three  exceptions, 
namely  (i)  that  the  special  parameters  are  now  only  ny  and  n2  since  surfaces 
are  two-dimensional,  (ii)  that  it  is  inappropriate  to  have  direct  input  of  metric 
terms  as  a  howstored  option,  and  (iii)  the  namelist  names  are  different. 


4.1.2  DE  and  SF 

Tags  DE  and  SF  have  no  special  parameters  and  allow  howstored  only  to  take 
the  value  file.  If  howstored  is  not  set,  then  default  values  are  set  for  variables 
in  the  namelist  that  follows  in  the  case  of  DE.  For  SF,  depending  on  class, 
parameters  are  set  that  otherwise  have  no  associated  tags. 
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4.2  List  of  Labels  Tags 

4.2.1  BL  and  PA 

Tags  BL  and  PA  define  a  multiblock  and  a  surface  patch  respectively.  The 
format  for  BL  is  specified  in  Section  3.2,  and  that  for  PA  differs  only  in  that 
the  label  following  the  patch  identifier  is  always  sameas,  thus  an  example  of  a 
record  is: 

PAupalysameasuperf ect.conductorl 

where  the  patch  physics  associated  with  label  perfecLconductorl  has  already  to 
have  been  defined  using  an  item  tagged  with  PP. 

4.2.2  OR 

The  tag  OR  is  followed  by  a  list  of  block  labels  in  the  order  of  assembly,  so  that 

e.g. 

ORublzublcybla 

indicates  that  information  associated  with  block  biz  supersedes  information  as¬ 
sociated  with  block  blc  which  in  turn  supersedes  information  about  block  bla. 
Equivalently  think  of  putting  biz  on  the  workbench,  then  attaching  blc  and  then 
bla. 

5  Class  Descriptions 

5.1  Block  Geometry 

5.1.1  Rectbk 

This  is  the  name  of  the  base  subroutine  used  to  mesh  orthogonal  geometries, 
that  have  either  a  Cartesian  reference  coordinate  system  (MCTYPE  =  1),  which 
is  the  default  case,  or  a  cylindrical  polar  reference  coordinate  system  (MCTYPE 
=  2).  Twelve  edge  curves  defined  by  discrete  points  along  their  length  are 
required  to  define  the  block.  MTYPE  is  zero  for  uniform  spacing  of  the  grid  in 
curvilinear  coordinates.  The  variables  in  NAMELIST  RECTBK  are 

MCTYPE  reference  coordinate  type 

MTYPE(3)  input  meshtype  selector 

MDEDGE  dimension  of  EDGES  array 

EDGES(MDEDGE)  array  of  edge  curves 

5.1.2  Genblk 

This  is  the  name  of  the  base  subroutine  used  to  mesh  nonorthogonal  geometries, 
that  have  either  a  Cartesian  reference  coordinate  system  (MCTYPE  =  1),  which 
is  the  default  case,  or  a  cylindrical  polar  reference  coordinate  system  (MCTYPE 
=  2).  Twelve  edge  curves  defined  by  discrete  points  along  their  length  are 
required  to  define  the  block.  MTYPE  is  zero  for  uniform  spacing  of  the  grid  in 
curvilinear  coordinates.  The  variables  in  NAMELIST  GENBLK  are 
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MCTYPE  reference  coordinate  type 

MTYPE(3)  input  meshtype  selector 

MDEDGE  dimension  of  EDGES  array 

EDGES(MDEDGE)  array  of  edge  curves 

5.1.3  Parallelepiped 

A  parallelepiped  is  specified  by  the  lengths  of  its  three  sides  and  the  (spherical) 
polar  angles  of  two  of  them,  the  third  side  being  assumed  to  be  along  the  polar 
axis  z.  The  variables  in  NAMELIST  PIPED  are 

RSIDl  length  of  parallelepiped  side  1 

RSID2  length  of  parallelepiped  side  2 

RSID3  length  of  parallelepiped  side  3 

THETl  polar  angle  of  paraUelepiped  side  1 

THET2  polar  angle  of  parallelepiped  side  2 

RPHIl  azimuthal  angle  of  paraUelepiped  side  1 

RPHI2  azimuthal  angle  of  paraUelepiped  side  2 

5.1.4  Polar  with  regular  meshing 

This  geometry  uses  a  cybndrical  polar  reference  coordinate  system.  The  vari¬ 
ables  in  NAMELIST  POLREG  are 


RADINR 

inner  radius 

RADOUT 

outer  radius 

THEMAX 

maximum  theta  (degrees) 

AXMAX 

maximum  along  cyl  polar-axis 

5.1.5  Variable  spaced  polar  lattice 

This  geometry  uses  a  cybndrical  polar  reference  coordinate  system.  The  vari¬ 
ables  in  NAMELIST  POLVAR  are 

RADINR  inner  radius 

RADOUT  outer  radius 

THEMAX  maximum  theta  (degrees) 

AXMAX  maximum  along  cyl  polar-axis 

5.1.6  Variable  spaced  cubic  lattice 

The  variables  in  NAMELIST  CUBVAR  are 

RXMAX  maximum  along  x-axis 

RYMAX  maximum  along  y-axis 

RZMAX  maximum  along  z-axis 
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5.1.7  Polar  mesh  segment 

The  variables  in  NAMELIST  POLSEG  are 


RADINR  inner  radius 

RADOUT  outer  radius 

THEMIN  minimum  theta  (degrees) 

THEMAX  maximum  theta  (degrees) 

AXMIN  minimum  along  cyl  polar-axis 

AX  MAX  maximum  along  cyl  polar- axis 


5.1.8  Polar  to  rectangular  transition 

The  variables  in  NAMELIST  POLRCT  are 

RADCUR  radius  of  curved  edge 
RADSTR  radius  of  straight  edge 

THEMIN  minimum  theta  (degrees) 

THEMAX  maximum  theta  (degrees) 

AXMIN  minimum  along  cyl  polar-axis 

AXMAX  maximum  along  cyl  polar-axis 


5.1.9  Regular  cubic  lattice 

The  variables  in  NAMELIST  CUBREG  are 

RXMAX  maximum  along  x-axis 

RYMAX  maximum  along  y-axis 

RZMAX  maximum  along  z-axis 


5.1.10  Quadrilateral  cylinder 

This  geometry  is  the  extrusion  of  a  quadrilateral  by  a  length  RZMAX.  In  a  plane 
z  =  constant,  the  corners  of  the  quadrilateral  are  in  clockwise  order:  (0,0), 
(x2,2/2)  and  {RY M AX,  RX M AX).  The  variables  in  NAMELIST 
QUADRI  are 


RXMAX 

RYMAX 

RZMAX 

XQUADl 

XQUAD2 

YQUADl 

YQUAD2 


maximum  along  x-axis 
maximum  along  y-axis 
maximum  along  z-axis 
xl  point  for  quadrilateral  cylinder 
x2  point  for  quadrilateral  cylinder 
yl  point  for  quadrilateral  cyhnder 
y2  point  for  quadrilateral  cylinder 
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5.2  Block  Physics 

5.2.1  Uniform 

The  variables  in  NAMELIST  UNIFRM  acre 

EPSR  relative  permittivity  in  uniform  block 

RMUR  relative  permeability  in  uniform  block 

5.3  Patch  Physics 

5.3.1  Perfect  conductor,  Applied  field,  Resistive  wall  and  Polar  axis 

These  classes  share  the  same  NAMELIST.  In  the  case  of  an  applied  electric 
displacement,  the  components  of  field  are  relative  to  the  coordinates  of  the 
block;  note  that  only  in  some  circumstances  is  d  constant.  STHETA  corresponds 
to  6  in  the  Annex  on  boundary  conditions.  The  variables  in  NAMELIST  BCS 
are 

SURFZ  surface  impedance 

DAPLYA(3)  applied  electric  field 

STHETA  lagging  factor  for  res.  wall  be 


5.4  Miscellaneous 

5.4.1  Free  parameters 

This  provides  the  opportunity  to  designate  a  specific  block  to  provide  the  ori¬ 
gin  of  coordinates  (using  those  variables  ending  in  GLB),  to  set  parameters  to 
control  a  run  (NRUN,  MDPART,  NOPSEL  and  NXPTDD),  to  set  initial  condi¬ 
tions,  and  to  control  diagnostic  output  as  described  in  the  Annex  on  3DP1C 
diagnostics.  DTMUL  is  the  Courant  number  Co  in  the  notation  of  the  Annex 
on  stability  and  dispersion.  MDPART  =  1  is  needed  to  prevent  the  introduc¬ 
tion  of  particles.  NINIT  corresponds  to  an  initial  uniform  magnetic  field,  an 
initial  cylindrical  TE  mode,  a  cylindrical  TM  mode,  a  rectangular  TE  mode, 
a  rectangular  TM  mode  or  a  magnetostatic  (PPM)  mode,  depending  whether 
it  takes  the  values  1,  2,  3,  4,  5  or  6.  The  field  components  are  set  by  BUNI 
and  mode  parameters  are  set  by  variables  ending  in  MODE.  Clearly  only  simple 
device  geometries  are  allowed,  and  the  3-coordinate  must  correspond  to  2.  The 
variables  in  NAMELIST  SF  which  have  not  already  been  described  in  Ref[4] 
are 


NPAR(MAXPAR) 

RPAR(MAXPAR) 

CBLGLB*32 

RFIGLB 

THEGLB 

XYZGLB(3) 


integer  parameter  defaults  on  tagged  line 
real  parameter  defaults  on  tagged  line 
block  with  designated  origin  and  rotation 
spherical  phi  rotn  of  named  block 
spherical  theta  rotn  of  named  block 
physical  position  of  named  block 
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DTMUL 

timestep  multiplier 

NRUN 

number  of  time  steps 

MDPART 

Dimension  of  PARTicle  coordinate  arrays 

NOPSEL 

Select  Output  Sequence 

NXPTDD 

select  0.4  EXPERT  diagnostic  dump 

AM0DE(3) 

box  size  for  modes 

BUNI(3) 

uniform  initial  B 

CM0DE(2) 

mode  amplitudes 

FMODE 

frequency  of  mode 

ZMODE 

impedance  of  mode 

NBEXT 

use  initial  B  as  external  B 

NINIT 

initial  field  type 

NM0DE(3) 

mode  numbers 
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Abstract 

This  paper  describes  three  dimensional  body-fitted  Particle- In-Cell  (PIC)  soft¬ 
ware  for  modelling  the  time  evolution  of  interactions  between  electromagnetic  (EM) 
fields  and  flows  of  relativistic  charged  particles.  A  description  is  given  of  the  phys¬ 
ical  model,  the  numerical  scheme  and  the  software.  The  performance  achieved  by 
the  software  on  parallel  computer  architectures  demonstrates  the  potential  of  this 
code  for  large  scale  time  domain  electromagnetic  and  electrodynamic  calculations. 

1  INTRODUCTION 

The  software  described  in  this  paper  offers  powerful  new  capabilities  for  electromagnetic 
PIC  modelling.  It  was  designed  for  three  dimensional  modelling  of  microwave  tubes  and 
microwave  transmission  where  the  interaction  of  electromagnetic  waves  with  charged 
particle  flow  is  important  [1,  2],  although  it  can  equally  well  be  applied  to  other  time 
dependent  problems  normally  tackled  by  finite  difference,  time  domain  (FDTD)  codes. 
Volume  filling  lossy  dielectric,  magnetic  media,  conducting  grids  and  surfaces,  surface 
physics  models  and  external  circuit  coupling  can  all  be  incorporated  into  the  simulations. 

It  differs  from  established  PIC  codes  [5,  6,  7,  8]  for  microwave  modelling  primarily  in 
that 

1.  it  uses  a  body  fitted  finite  element  net  rather  than  orthogonal  finite  difference  nets, 

2.  the  finite  element  (FE)  derivation  gives  tensor  field  equations  which  carry  over 
the  simplicity  of  the  Yee  FDTD  algorithm  [9]  for  the  EM  calculation  to  general 
geometry  nets, 

3.  the  FE  derivation  gives  a  charge  conserving  current  assignment  scheme  on  the 
general  geometry  nets, 

4.  it  was  designed  ab  initio  to  use  a  multiblock  element  net  decomposition  to  use 
efficiently  message  passing  in  distributed-memory  parallel  computers. 

The  next  section  presents  the  physical  model,  followed  in  Section  3  by  an  outline 
of  the  numerical  method.  Major  problems  are  the  description  of  objects  and  compact 
data  storage  for  the  computations;  these  issues  are  addressed  in  the  Section  4.  Parallel 
computer  implementation  is  treated  in  Section  5 


2  PHYSICAL  MODEL 


The  basic  physics  involved  is  well  described  by  Maxwell’s  equations  for  the  electromag¬ 
netic  fields  and  Vlasov’s  equation  (possiblj'  relativistic)  for  the  dynamics  of  the  charged 
particles.  The  charged  particles  may  be  electrons,  ions,  or  both,  although  the  illustrations 
given  below  are  for  electrons  only.  The  physical  model  is  completed  by  adding  descrip¬ 
tions  of  the  emission  and  absorption  of  radiation  and  of  particles  from  boundaries.  In 
microwave  tubes  these  boundaries  may  have  complex  geometries  and  may  be  coupled  to 
external  circuits.  Within  the  devices  may  be  lossy  dielectric  and  magnetizable  materials. 

The  mathematical  model  is  the  relativistic  Maxwell-Vlasov  set.  If  the  distribution 
function  is  represented  by  a  set  of  sample  points  (i.e.  “superparticles”) 

/(x,p,t)  =  Y^N,S{x-Xi{t))6{p  -  pi{t))  (1) 

t 

where  (x,,  p,);  i  —  [l,Ay  are  the  coordinates  of  the  Np  superparticles,  each  of  mass  M 
=  NsItiq  and  charge  Q  =  Ngq,  then  the  Maxwell-Vlasov  set  may  be  written  in  terms  of 
the  action  integral 

,  f  ,  ,  fD  E-B  H  ^  .A  f^Mc\ 

I^Jcltdry - - - p(l)  +  y  Aj  -  J  (2) 

where  B  =  V  x  A  and  E  =  —  V^  —  A.  Treating  I  as  a  functional  of  the  vector  potential 
A,  the  scalar  potential  4>  and  particle  coordinates  {x,  }  leads  to  Euler-Lagrange  equations 
which  reduce  to  Maxwell’s  equations  and  the  relativistic  equations  of  motionflO]. 

2.1  Lossy  dielectrics 

General  dielectric  and  magnetic  media  may  be  added  to  the  model  by  introducing  polar¬ 
isation,  P,  and  magnetisation,  M,  terms  into  the  action  integral 

/  =  ...  +  J  didriM  •  B  +  P  •  E)  (3) 

and  substituting  for  M  and  P  from  the  resulting  evolutionary  equation  using  their  con¬ 
stitutive  relationships.  The  non-lossy  part  gives  the  dielectric  function  relating  D  to  E, 
and  the  lossy  parts  give  the  magnetisation  and  conduction  currents,  j/:- 

D  =  cE,  j(  =  —V  X  (3B  +  crE  (4) 

It  has  been  found  advantageous  in  simulations  to  use  a  small  magnetisation  current  to 
control  numerical  noise.  This  allows  simulations  to  be  performed  with  fewer  superparti¬ 
cles  than  would  otherwise  be  possible. 

2.2  External  Circuits 

The  full  wave  description  of  the  microwave  tube  is  linked  to  a  lumped  circuit  approxi¬ 
mation  of  power  supplies  and  extraction  waveguides  by  coupling  elements  at  the  surfaces 
of  the  active  computational  domain  to  external  circuit  elements.  The  external  circuit 


elements  have  no  ph3^sical  dimension,  but  they  provide  a  relationship  between  surface 
current,  /,  and  tangential  electric  fields,  E.  Thus,  a  purely  resistive  circuit  element 
would  give  E  =  ZI  and  a  circuit  with  resistance  and  inductance  (or  capacitance)  would 
give 

where  coefficients  a,  (3,^,6  are  chosen  to  describe  the  particular  external  circuit.  Simi¬ 
larly,  a  second  order  differential  equation  is  used  for  coupling  to  an  L-R-C  circuit  and  so 
forth. 

2.3  Electron  Emission 

Electron  emission  from  boundary  surfaces  is  treated  using  standard  Monte-Carlo  proce¬ 
dures.  Beam  injection  selects  superparticle  positions  and  momenta  to  fit  the  incoming 
beam  distribution  function.  Space  charge  limited  emission  at  cathodes  is  modelled  by 
introducing  sufficient  free  surface  space  charge  (in  the  form  of  superparticles)  to  bring 
the  normal  field  at  the  surface  to  zero.  Secondarj'  electrons  at  boundary  surfaces  are 
generated  and  emitted  to  fit  the  chosen  energj'  and  direction  distribution  in  response  to 
the  locations  and  momenta  of  impacting  superparticles. 


3  NUMERICAL  METHOD 


The  discrete  equations  are  obtained  using  the  virtual  particle  method  [4],  with  fields  in 
the  general  geometry  version  being  represented  by  their  tensor  components.  These  com¬ 
ponents  are  approximated  by  finite  elements  in  both  space  and  time,  and  the  distribution 
functions  are  represented  by  sets  of  sample  points  (‘superparticles’).  The  finite  element 
and  sample  approximations  are  substituted  into  the  action  integral,  whence  variations 
with  respect  to  potentials  lead  to  discrete  forms  of  Maxwell’s  equations  and  prescriptions 
for  charge  and  current  assignment.  Current  is  assigned  from  ‘virtual  particles’  placed 
at  points  specially  interpolated  between  positions  at  successive  time  levels,  a  procedure 
which  automatically  leads  to  charge  conservation.  In  all  cases,  we  have  assumed  the 
lowest  order  conforming  elements. 

In  general  curvilinear  coordinates  {x^,  x^)  the  Action  integral  is 

1  =  jdt  dx^  dx^  dx^  |i(E,-  d'  -  H,  V)  +  P  Ai  -  Q(^|  (6) 


where 


d(f>  dAk  •  iikdAk 

ail 


d'  =  toy/gg'^Ej,Hi  = 


1 


9ii^ 


The  ‘superparticle’  current  is  given  bj”^  the  sum  over  all  particles: 

=  S  %d{x:^  -  d:l)S{x^  -  il)6{x^  -  xl)x 


(7) 

(8) 

(9) 


3,1  Field  Equations 


Treating  Action  I  as  a  function  of  the  finite  element  approximations  for  (f)  and  A,,  then 
taking  variations  with  respect  to  the  nodal  amplitudes  gives  the  element  contributions  to 
the  discrete  approximations  to  Maxwell’s  equations.  In  general,  the  form  of  the  assembled 
equations  depends  on  the  number  of  elements  and  boundaries  adjacent  to  the  node  in 
question.  However,  for  an  internal  node  where  the  element  net  is  topologically  equivalent 
to  a  cubic  lattice,  the  assembled  equations  take  the  familiar  form: 


d,W  =  0 

(10) 

-  I’, 

,  d,6'  =  Q 

(11) 

where  the  symbol  d  denotes  a  centred  difference.  These  equations  in  tensor  quantities 
b',d’,  E/.,H^.,I‘  and  Q  take  the  same  form  as  the  Yee  FDTD  scheme  [9]  for  cartesian  field 
components,  but  are  ideniicol  in  any  coordinate  systeni.  This  provides  great  simplifica¬ 
tion. 

3.2  Constitutive  Equations 

Geometrical  and  material  (permeability  and  permittivity)  informat  ion  appears  only  in  the 
constitutive  relations.  The  weak  approximations  to  Eqs  (8),  with  lumped  mass  matrices 
give  the  simplest  explicit  expressions  for  E,-  and  H;  of  the  form: 

E,=GSd>;  H,=Gjb^  (12) 

Elements  of  the  symmetric  tensors  Gfj  and  Gfj  are  sparse  matrices.  More  general  per¬ 
mittivities  and  permeabilities  are  handled  by  replacing  Cq  and  //q  by  scalar  or  tensor 
functions  in  the  computation  of  and  Gfj,  respectively. 

3.3  Assignment 

The  formulae  for  current  and  charge  assignment  arising  from  the  variational  formulation 
are  also  coordinate  independent.  Let  the  finite  element  test  function  approximations  to 
the  potentials  be  d>  =  "  and  A,  =  A(,)1T(.),  where  and  A  are  nodal  amplitudes,  and 

sums  are  implied  over  element  nodes  (cf  Ref  [4]).  Then  variations  of  Eq(6)  yield: 

charge  assignment  :  Q  =  f  dtJ2qpU{{xl,  xl)it),t)  (13) 

V 

current  assignment  :  /'  =  /  ^  9p.Tj,‘^iy(,)((.fJ,  xj,  .Tp)(t),  t)  (14) 

•'  V 

In  two  dimensions,  the  integrals  for  Q  and  /'  are  evaluated  in  exactly  the  same 
manner  as  described  in  Ref  [4].  The  only  difference  is  that  the  cartesian  coordinates 
(a:,?/)  are  replaced  by  the  curvilinear  (a^^,.f^).  This  generalises  straightforward!}"  to  three 
dimensions,  but  then  requires  assignment  from  two  rather  than  one  ‘virtual  particle’  per 
segment  per  element  to  ensure  exact  charge  conservation. 


3.4  Charge  Conservation 

Charge  and  current  assignment  are  linear  operations,  so  bj^  linear  superposition,  conser¬ 
vation  for  one  trajectory  moving  through  a  single  time  step  implies  the  same  for  the  sum 
of  all  trajectories.  Summing  all  contributions  to  currents  from  a  single  particle,  gives 
the  same  as  the  difference  of  charge  at  start  and  end  times,  ie  the  linear  quadrilateral 
(hexahedral  in  three  dimensions)  element  case  satisfies 

dtQ  =  -dkl^  (15) 

where  the  symbol  d  denotes  a  centred  difference  arising  from  assembling  the  finite  element 
contributions,  and  Q  and  are  nodal  amplitudes.  Equations  of  the  form  (15)  can  be 
shown  to  be  generallj^  satisfied  for  VP  algorithms. 


3.5  Particle  Dynamics 

Particle  equations  of  motion  are  integrated  bj'^  first  updating  the  curvilinear  momentum 
components,  measured  at  the  current  particle  positions  using 

-  pV'"  =  +  n'-"+'/2)6V2  (16) 

where  are  contravariant  velocity  components  and  tkim  is  the  permutation  symbol.  The 
Christoffel  symbol  term  does  not  appear  as  Eq(16)  is  evaluated  at  a  fixed  position. 

Positions  are  updated  to  new  timelevel  +  1  and  the  momentum  components  at  the 
new  positions  are  computed  using  a  predictor-corrector  scheme  on 


and 

p'’n+i/2(a'”+i)  =  e'(.r"+^)  •  efc(.T")p''’”+'/^(a:")  (18) 

In  cartesians,  this  scheme  reduces  to  the  usual  Lorentz  force  leapfrog  integration  scheme 

[11]. 


4  SOFTWARE  STRUCTURE 

The  problem  in  designing  the  software  is  to  devise  a  scheme  where  complex  shaped  objects 
can  be  simulated  efficiently  on  a  distributed  memory  parallel  computer.  The  solution  we 
have  adopted  is  to  divide  the  object  being  modelled  into  a  set  of  curvilinear  hexahedral 
blocks,  where  each  block  is  equivalent  to  a  cube.  Each  block  has  its  own  fields  and 
particles  and  communicates  with  other  blocks  and  boundaries  only  via  its  patch  buffers. 

Global  data  describes  the  location  and  logical  connection  of  the  blocks,  and  is  common 
to  all  the  processors  on  which  the  program  resides.  Local  data  describes  fields  and  particle 
coordinates,  and  is  different  in  each  processor’s  copy  of  the  program.  If  blocks  reside  on 
the  same  processor,  their  patch  data  is  exchanged  by  memory  to  memory  copies,  but  if 
they  are  on  different  processors,  message  passing  is  used. 

Data  initialisation  is  performed  as  follows:  A  device  is  specified  by  compact  descrip¬ 
tions  of  the  building  blocks,  their  connectivity,  and  the  boundary  condition  patches  on 


surfaces.  Input  validation  checks  that  the  surface  of  every  block  is  completely  covered 
by  patches  and  that  the  connections  are  physically  realisable.  Restricting  boundary  con¬ 
ditions  to  block  surfaces  greatlj'  simplifies  their  application,  but  does  not  constrain  load 
balancing.  Several  blocks  may  be  assigned  to  one  processor,  or  larger  blocks  may  be 
subdivided  to  balance  loading  across  processors. 

The  volume  within  each  block  is  subdivided  into  a  lattice  of  elements  using  transfinite 
interpolation  [12].  In  the  curved  space  in  which  the  equations  are  solved,  the  lattice 
becomes  a  rectangular  brick  of  cubical  elements,  leading  to  simple  addressing  and  fast 
code.  We  use  a  data  compression  scheme  which  makes  the  storage  needed  for  basis 
vectors  and  metrics  negligible  for  most  blocks. 

5  PARALLEL  PERFORMANCE 


1  10  100  1000 
Number  of  Processors,  p 

Figure  1:  Temporal  performance  of  the  MIMD  PIC  code  for  three  problem  sizes  as  a 
function  of  the  number  of  processors  used,  on  a  16  processor  Intel  iPSC/860.  Solid 
symbols  are  for  per-patch  code,  open  symbols  for  per-processor  code.  The  dotted  lines 
are  the  theoretical  ideal  behaviour  scaled  from  the  one-processor  performance  for  the 
8-block  case. 

The  message  passing  model  used  in  our  program  matches  the  structure  of  the  new  gener¬ 
ation  of  massively  parallel  processors  (MPPs)  such  as  the  Intel  Paragon  and  Meiko  CS-2, 


and  is  also  implemented  efficiently  on  the  traditional  parallel  vector  shared-memory  com¬ 
puters  such  as  the  Cray  Y-MP  and  C-90.  Our  code  currently  uses  the  Intel  NX/2  library 
of  communication  subroutines  (mostly  CSEND  and  CRECV)  because  these  are  imple¬ 
mented  with  minimum  overhead  on  many  computers  and  are  widely  used.  Conversion 
of  the  code  to  PVM  [13]  or  the  emerging  MPI  standard  [14]  should  be  straight  forward 
since  message-passing  calls  are  confined  to  a  few  subroutines.  We  did  not  use  PVM  for 
our  implementation  because  of  the  very  high  startup  time  for  message  send  on  most 
implementations  [15]. 

Both  two  and  three  dimensional  test  cases  have  been  used  to  evaluate  parallel  perfor¬ 
mance.  Shown  in  figure  1  are  measurements  for  three  sizes  of  three  dimensional  electron 
plasma  calculations  run  using  up  to  64  processors  on  an  Intel  iPSC/860.  Cubical  plasma 
computations  comprising  8,  64  and  512  blocks  of  64  elements  were  used.  Each  block  was 
initialised  with  512  particles.  For  each  problem  size  the  performance  is  measured  for  a 
variety  of  numbers  of  processors  in  order  to  test  the  scaling  behaviour.  We  measure  the 
performance  using  an  absolute  performance  metric,  namely  the  Temporal  performance 
[17]  expressed  in  units  of  timestep  per  second  (tstep/s),  because  this  is  the  figure  of  merit 
that  the  user  of  any  simulation  program  wishes  to  maximise. 

Two  performance  curves  are  given:  the  solid  symbols  show  the  results  for  the  per- 
patch  code  which  sends  a  separate  message  for  each  glue-patch,  whilst  the  corresponding 
open  symbol  shows  the  better  performance  that  is  obtained  if  the  patch  messages  are 
combined  so  as  to  send  only  a  single  multi-patch  message  between  each  pair  of  connected 
processors,  the  per-processor  code.  The  smallest  8-block  problem  can  only  use  8  pro¬ 
cessors  at  the  most  (one  block  per  processor),  and  all  possible  points  are  measured  and 
shown.  The  increase  of  performance  with  number  of  processors  is  very  close  to  the  ideal 
linear  rise  proportional  to  the  number  of  processors,  which  is  shown  by  the  dotted  line 
rising  at  45  degrees.  Ideally,  not  only  do  we  aim  to  achieve  ’p’  times  the  performance 
of  a  single  processor  on  a  given  problem  when  using  p-processors,  we  also  hope  to  be 
able  to  solve  a  problem  that  is  ’p’  times  bigger  in  the  same  time.  The  other  dotted  ideal 
performance  lines  express  this  relation,  and  are  scaled  from  the  one-processor  point  for 
N't,  =  8-  The  measured  results  for  the  per-processor  code  for  the  other  two  problem  sizes 
follow  this  scaling  relationship  almost  exactly.  In  fact,  the  worst  Efficiency  (percentage 
deviation  from  the  ideal  dotted  lines)  for  the  per-processor  code  is  76%,  which  is  seen 
for  the  64-block  case  on  64  processors.  We  may  conclude  therefore  that  the  scaling  per¬ 
formance  of  the  code  is  good  up  to  64  processors.  We  anticipate  similar  scaling  to  be 
maintained  for  larger  problems  and  the  larger  number  of  processors  available  on  the  Intel 
Paragon. 
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ABSTRACT 

A  relativistic  electromagnetic  PIC  algorithm  for  general  three  dimensional  ge¬ 
ometries  is  described.  The  correct  combination  of  co-  and  contravariant  field 
components  and  of  weighting  yields  simple  coordinate  invariant  numerical  pre¬ 
scriptions.  The  use  of  isoparametric  hexahedral  elements,  generated  by  trans- 
finite  interpolation,  and  multiblock  decomposition  leads  to  algorithms  ideally 
suited  to  MIMD  computers,  as  demonstrated  by  the  speedup  obtained  on  an 
iPSC/860. 

1  Introduction 

Virtual  Particle  (VP)  particle-mesh  algorithms  are  now  established  as  an  effective  approach  to 
obtaining  numerical  schemes  for  solving  the  relativistic  Maxwell- Vlasov  equations  [1,  2,  3,  4]. 
Unlike  conventional  Particle-in-Cell  (PIC)  schemes,  they  are  derived  using  finite  elements  in 
both  space  and  time.  Current  is  assigned  from  ‘virtual  particles’  placed  at  points  specially 
interpolated  between  positions  at  successive  time  levels,  a  procedure  which  automatically 
leads  to  charge  conservation.  Existing  VP  implementations  use  rectangular  finite  elements 
in  two  dimensional  Cartesian  and  polar  geometries.  Only  a  restricted  class  of  device  is  well 
modelled  in  such  circumstances,  leading  to  the  need  to  implement  VP  in  more  complex 
geometries.  Section  2  of  this  paper  shows  how  the  VP  method  extends  to  general  three 
dimensional  body-fitted  elements  [5],  whilst  Section  3  explains  how  blocks  of  these  elements 
are  assembled  to  model  realistic  devices,  and  shows  the  speedup  achievable  on  distributed 
Memory  MIMD  architecture  computers. 


2  Variational  Formulation 


In  general  curvilinear  coordinates  (x^, 

I  =  J  dt  dx^  dx^  dx^^/g 


x^,  x^)  the  action  integral  may  be  written 

{^{Ei  -  Hi  B^)  -f  r  Ai  -  +  Ik 


where  Ik  is  the  kinetic  Lagrangian  and 

d(f)  dAk 
~  dx>^  dt 


(1) 

(2) 
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_  t”'  dAt 

y/g  dxJ 

(3) 

D'  =  =  e^g'^Ej 

(4) 

//.  =  —  =  -g.jB^ 

Po  Po 

(5) 

If  we  further  assume  that  the  distribution  function  is  represented  by  a  set  of  sample  points 
(ie  ‘superparticles’),  then  the  source  terms  in  the  field  Lagrangian  become 


T  \/9 

and  the  kinetic  Lagrangian  term  becomes 


(6) 

(7) 

(8) 


The  sums  in  p  are  over  particles,  each  with  charge  qp.  The  metric  tensor  elements  gij  can  be 
computed  from  the  relationship  between  Cartesian  and  the  general  curvilinear  coordinates 
i.e.  given 

x'  =  x^)  (9) 


then 


dx^ 


(10) 


where  g  =  ||<7,j|l  and  gijg^^  —  Treating  I  as  a  functional  of  the  vector  potential  Ai,  the 
scalar  potential  (j)  and  particle  co-ordinates  {Xpj  leads  to  Euler-Lagrange  equations  repre¬ 
senting  Maxwell’s  ecjuations  and  relativistic  particle  motion. 

Substituting  test  function  approximations  for  (j),  Ai  and  Xp  and  varying  with  respect  to 
the  nodal  amplitudes  yields  a  finite  element  approximation  to  the  Maxwell-Vlasov  equations. 
Now,  if  we  introduce  b‘  =  y/gB\  d'  =  y/gD\  J'  —  y/gj'  and  Q  =  y/gp,  Maxwell’s  equations 
become 


dt  dxi'  dx' 


dt  dx^  '  dx' 


Q. 


(11) 

(12) 


where  e'-^^  is  the  permutation  symbol.  This  suggests,  and  our  analysis  confirms,  that  the  dis¬ 
crete  VP  equations  in  the  quantities  6',  d',  i/,,  J'  and  Q  are  identical  in  any  co-ordinate 

system.  Geometrical  information  appears,  along  with  the  permeability  and  permittivity  ten¬ 
sors,  only  in  the  constitutive  relations  relating  H  to  B  and  D  to  E. 
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3  Mesh  Generation  and  MIMD 


A  general  curvilinear  coordinate  system  arises  naturally  when  we  seek  to  represent  awkwardly 
shaped  devices.  Transfinite  interpolation  provides  a  coordinate  system  for  any  simply  con¬ 
nected  surface  (in  2-D)  or  volume  (in  3-D)  provided  the  bounding  curves  or  surfaces  are  not 
too  convex  or  concave.  Realistic  microwave  sources  involve  cavities  and  waveguides  of  cylin¬ 
drical  and  rectangular  cross-section.  Separate  subdomains  or  blocks  are  therefore  introduced 
to  handle  the  various  components  and  also  axial  effects. 

The  multiblock  subdivision  of  the  computational  domain  minimises  the  amount  of  global 
data  and  interprocessor  message  passing,  and  simplifies  load  balancing  across  processors. 
Each  slave  block  only  requires  data  from  its  neighbours,  and  the  master  control  program 
only  requires  information  about  the  block  surfaces  (‘glue  patches’)  which  join  the  slave  blocks 
together.  This  arrangement  offers  the  prospects  of  large  computational  intensitj'  and  a  weak 
Amdahl  limit  to  speedup  on  distributed  memory  MIMD  computers.  Moreover,  the  simple 
logical  cube  addressing  within  each  block  leads  to  fast  serial  processing. 

We  can  further  demand  that  boundary  conditions  only  apply  at  the  surfaces  of  blocks;  this 
eliminates  the  addressing  problems  in  embedding  surfaces  within  blocks,  and  allows  surface 
data  to  be  passed  to  the  control  program  through  the  ‘glue  patch’  tables.  Further  saving  of 
computer  storage  and  time  arise  from  keeping  metric  information  and  material  property  data 
only  in  those  blocks  where  they  are  needed.  When  many  small  blocks  are  used  to  describe  a 
complex  object,  load  balancing  is  achieved  by  assigning  several  blocks  to  one  processor. 

3.1  Test  Problem 

Figure  1  shows  results  obtained  using  an  iPSC/860  hypercube  for  a  small  test  problem 
involving  a  planar  MILO  configuration.  The  device  is  modelled  using  20  multiblocks  that 
contain  either  192  or  128  elements;  the  total  number  of  particles  emplo5'’ed  is  under  5000.  The 
run  is  chosen  to  reflect  the  general  characteristics  of  a  production  calculation,  except  that  it 
is  terminated  after  only  a  few  (100)  timesteps.  It  also  differs  from  the  normal  MILO  runs  in 
that  the  device  is  initially  filled  with  a  nonzero  electric  and  magnetic  field.  The  outcome  is 
that  electrons  are  emitted  from  the  whole  length  of  the  cathode,  and  there  is  a  strong  initial 
transient  where  the  electrons  fill  only  part  of  the  device  volume. 

The  results  of  this  and  similar  benchmark  computations  show  encouraging  speedup,  even 
for  modestly  sized  calculations.  Over  80%  ideal  linear  speedup  is  achieved  even  for  very 
non-uniform  electron  distributions,  cf.  Figure  1.  The  present  implementation  bases  message 
passing  on  a  patch  to  patch  basis.  On  machines  with  high  interprocessor  message  passing 
latency,  further  speedup  would  result  from  presorting  the  patches  and  performing  message 
passing  on  a  processor  by  processor  basis. 


4  Final  Remarks 

Our  preliminary  VP  benchmark  computations  support  the  contention  that  the  multiblock 
decomposition  allowed  by  these  algorithms  is  ideal  for  distributed  memory  MIMD  computers. 
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Number  of  Processors 

Figure  1:  Temporal  Performance  in  timesteps  per  second  (tstep/s) 

The  electromagnetic  calculation  can  be  statically  load  balanced,  has  little  message  passing 
between  blocks  and  small  scalar  overheads. 

We  are  currently  developing  general  geometry  software  based  on  the  scheme  described 
herein  in  order  to  extend  our  capability  for  microwave  tube  modelling.  However,  the  elegant 
body-fitting  electromagnetic  solver  part  of  our  PIC  scheme  can  equally  well  be  applied  to 
the  computationally  simpler  problems  where  free  charges  are  absent;  for  example,  to  study 
waveguides  and  exterior  scattering  problems. 
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